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modified  to  simplify  calculations  by  excluding  from  consideration 
those  minor  constituents  whose  total  contribution  to  the  mass  density 
is  less  than  1%."\ 


> Software -haSj^been  developed  to  support  studies  of  atmospheric  density 
variations  based  on  data  collected  by  satellite-borne  accelerometers.  , 

i  t 

Modifications  4wiveA been  made  to  several  standard  satellite  ephemeris 
codes.  The  changes  include:  (1)  adaptation  of  the  imbedded  atmospheric 
density  model  to  extend  its  validity  to  cases  of  extremely  high 
solar  activity;  and  (2),  improvement  in  the  accuracy  of  conversion 
between  mean  Keplerian  elements  and  position/  velocity  vectors., 

A  spurious  signal  originating  in  the  receiver  is  found  to  contaminate 
recorded  data  from  scintillating  trans-ionospheric  radio^waves  transmitted 
from  a  satellite  beacon.  Special  pre-processing  is  shown  to  eliminate 
the  effect  of  the  contaminant  on  statistical  analyses  of  the  scintillation 
data. — 

" An  HF  ducted  mode  of  propagation  is  to  be  supported  using  artificial 
ionospheric  irregularities.  One  terminal  is  to  be  ground-based;  the 
other,  satellite-borne.  A  standard  satellite  ephemeris  program  has 
been  augmented  to  provifS* time  histories  of  range  and  range-rate  for 
both  the  ducted  mode  anti 'the  classical  ionospheric  reflection  modes. 
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Atmospheric  Density  Models 


This  section  describes  recent  analysis  and  programming  performed  by 
Logicon,  Inc.,  to  develop  improved  upper  atmospheric  neutral  density 
models.  Knowledge  of  upper  atmospheric  neutral  densities  is  an 
important  requirement  for  the  understanding  and  analysis  of  many 
phenomena  under  study  at  AFGL.  Researchers  in  many  areas  desire  best 
estimates  of  composition  and  temperatures  as  inputs  to  models  and  data 
analysis  programs.  Prominent  examples  are  analyses  of  auroral, 
airglow,  and  ionospheric  measurements. 

The  continuing  need  for  accuracy  in  satellite  tracking  and  ephemerides 
prediction  also  results  in  the  need  for  improved  modeling  of 
thermospheric  density  variations.  Poorly  modeled  magnetic  storm  and 
local  time  variations  continue  to  limit  satellite  tracking  and 
prediction  accuracy,  particularly  at  lower  altitudes.  Operational 
models  used  in  this  application  must  be  not  only  accurate  but  also 
efficient.  The  computation  of  satellite  drag  requires  only  the  total 
mass  density,  as  opposed  to  the  composition.  Thus  operational  models 
should  be  formulated  to  compute  the  total  mass  density  directly  if 
possible,  or  to  limit  the  composition  to  the  most  significant 
components. 


Extensive  research  in  atmospheric  density  modeling  has  continued  at 
AFGL.  Theoretical  studies  of  local  time  variations  have  been 
conducted. ^  Section  1.1  describes  Logicon's  efforts  to  incorporate  the 
results  of  these  studies  into  easily  computable  formulations  of 
temperature  and  density  local  time  variations.  Section  1.2  discusses 
modifications  made  to  AFGL's  CADNIP/BADMEP^  system  for  satellite 
orbital  decay  analysis  and  prediction.  Some  of  these  modifications  are 
temporary  changes  to  an  existing  density  model  for  studies  at  extremely 
high  solar  activity.  In  addition  permanent  mod  if ict ions  have  been  made 
to  correct  an  inaccuracy  in  the  conversion  between  mean  elements  and 
position-velocity  and  to  permit  direct  input/output  of 
position-velocity.  Section  1.3  describes  8  modification  to  the 


Jacchia-Bass^  density  model  to  exclude  components  whose  total 
contribution  to  the  total  mass  density  is  less  than  IS.  Finally, 
Section  1.4  describes  software  developed  for  studies  of  atmospheric 
density  variations  using  satellite  accelerometer  data  collected  by  AFGL 
scient ists^*^»^. 

1.1  Atmospheric  Tides 

While  the  response  of  the  upper  atmosphere  to  geomagnetic  activity  and 
unmodeled  density  waves  continue  to  be  the  major  sources  of  error  in 
practical  neutral  thermospheric  models,  variations  with  local  time  are 
also  poorly  represented  in  most  models,  particularly  at  low  altitudes. 
This  is  because  the  behavior  switches  from  primarily  diurnal  (24  hour 
period)  at  high  altitude  to  semidiurnal  (12  hour  period)  at  low 
altitude, ^  while  much  of  the  data  used  to  constuct  the  models  are  at 
high  altitudes.  An  exception  to  this  is  the  MSIS  model, ^>9  which, 
however,  is  based  largely  on  midlatitude  data^.  Furthermore,  the  MSIS 
model  does  not  represent  well  the  atmospheric  response  to  geomagnetic 
activity,  since  it  uses  only  daily  averages  of  the  geomagnetic  activity 
index  rather  than  3-hour  averages. 

Work  by  Forbes  and  Garrett^  on  atmospheric  tides  encompasses  both 
theoretical  considerations  and  data,  in  that  certain  parameters  of  the 
model  are  calibrated  from  data.  In  particular,  the  semidiurnal  tide  in 
the  thermosphere  is  composed  of  two  parts:^  that  due  to  direct 
(in-situ)  excitation  by  solar  radiation  and  ion-neutral  momentum 
coupling,  and  that  due  to  upward  propagation  of  modes  excited  in  the 
mesosphere.  The  in-situ  portion  is  computed  by  solving  a  set 
linearized  equations  expressing  momentum  and  energy  conservation, 
continuity,  and  the  perfect  gas  law.  In  these  equations  the  input 
solar  energy  is  determined  by  requiring  the  exospheric  diurnal  tide, 
which  is  assumed  to  be  in-situ,  to  agree  with  data.  Fourier 
decomposition  in  local  time  of  the  heat  source  determines  the 
semidiurnal/diurnal  excitation  ratio,  which  together  with  the  diurnal 
excitation  determined  as  indicated,  provides  the  semidiurnal  heat 
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input.  The  contribution  due  to  the  upward  propagation  of  tidal  modes 
from  the  mesosphere  is  determined  by  subtracting  the  computed  in-situ 
tide  from  available  wind  and  temperature  measurements  and  fitting  the 
residuals  to  a  linear  combination  of  such  modes H. 

In  Section  1.1.1  will  be  discusssed  the  development  of  a 
computationally  practical  representation  of  the  total  density  annually 
averaged  semidiurnal  tide  based  on  the  Forbes-Garrett  model  just 
discussed.  This  requires  first  a  computation  of  the  tidal  variations  in 
composition,^  the  results  of  which  are  combined  with  a  suitable 
background  composition  model  to  compute  total  density  tides.  The 
resulting  semidiurnal  tide  is  then  fit  to  simple  functional 
representations  in  height  and  latitude  for  incorporation  into  existing 
semi-empirical  codes.  Simulation  of  diurnal  phase  and  amplitude 
variations^  is  also  discussed.  Efforts  to  develop  a  temperature  tidal 
representation  for  use  in  composition  models  such  as  the  MSIS  and 
3acchia-Bass  models  are  discussed  in  Section  1.1.2. 

1.1.1  Total  Density  Tides 

The  development  of  a  total  mass  density  tidal  model  begins  with  the 
computation  of  compositional  tides.  The  computation  was  carried  out 
following  Forbes'  formulation^*^  in  which  the  heiqht  variation  of  the 
tidal  perturbation  in  the  log  of  the  number  density  of  the  1th 
const itutent  is  governed  by  the  equations: 


{7*?  +  -Y1Bi01 
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,  *  ,  m 


where 


Ri  =  vinNi 


h  = 

Ni  = 

a1  = 
To  = 

T'  = 
Hio  = 
Hi*  = 
Di  = 

*io  = 
n  = 


v.V  = 
W  = 
-KN*  = 

j  = 


altitude 

number  density  of  ith  const itutent 

thermal  diffusion  coefficient  of  ith  constituent 

time  averaged  temperature 

tidal  perturbation  in  temperature 

time  averaged  scale  heiqht  of  ith  constituent 

tidal  perturbation  in  scale  height  of  ith  consitutent 

diffusion  coefficient  of  ith  constituent  through  major 

gas 

time  averaged  vertical  flux  of  ith  constituent 

zonal  wave  number  (1  for  diurnal  tide,  2  for  semidiurnal 

tide) 

rotation  rate  of  earth 
divergence  of  horizontal  velocity 
vertical  velocity  of  major  gas 

loss  rate  of  N-j  due  to  chemical  reaction  with  some 
constituent 

VT- 


Program  TIDEVAR  (Figure  1)  was  constructed  to  carry  out  these  ’• 

computations.  This  program  computes  the  tidal  perturbations  for  the 
constituents  O2,  0,  N2,  N,  He,  Ar,  and  H.  Atomic  nitrogen  has  been  ( 

included  for  possible  use  in  future  composition  studies,  although  it  5 
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PROGRAM  TIOVAR( InPUT , OUTPUT , TAPE  1 , T APE2 , TAPE3*520 , TAPE4-520 ) 
THERMOSPHERIC  COMPOSITION  TIDAL  VARIATION  PROGRAM 
FORBES, J.  (JOURN.  GEOPHYS.  RES..VOL.  83, NO.  A8 , 

PP.  3691-3698, 197B) 

DIFFUSIVE  EQUILIBRIUM  MAY  BE  ASSUMED, I.  E..  EOUATION  (3) 

OF  FORBES  IS  SOLVED  WITH  WI'-W  «  0. 

IF  DIFFUSIVE  EQUILIBRIUM  IS  NOT  ASSUMED, 

EOUATION  (8)  IS  SOLVED,  BUT  THE  CHEMISTRY  TERM  KNJ 1 
IS  NEGLECTED. 

INPUTS 

CARD  1  (LIST-OIRECTED) 

iOPT  -  BOUNDARY  CONDITION  INDICATOR  (SEE  X  ANO  Y  BELOW) 

HO  -  BOTTOM  REDUCED  HEIGHT(IN  SCALE  HEIGHTS  ABOVE  GROUND) 

HN  -  TOP  REDUCED  HElGHTl SCALE  HEIGHTS  ABOVE  GROUND) 

DH  -  REDUCED  HEIGHT  STEP  SIZE(SCALE  HEIGHTS) 

X  (COMPLEX).  Y 

FOR  IOPT  a  1  XsBOTTCM  FUNCTION  VALUE. 

Y  IGNORED 

FOR  IO?T«2  X a FUNCTION  VALUE 

AT  REDUCED  HEIGHT  Y 
THE  OTHER  BOUNDARY  CONDITION , 

THE  FIRST  DERIVATIVE  AT  THE  TCP, 

IS  SET  AUTOMATICALLY  FROM 
DIFFUSIVE  EQUILIBRIUM 
CARD  2  (LIST-DIRECTED) 

IYR  -  YEAR 
IMO  -  MONTH 
IDA  -  DAY 

XLAT  -  LATITUDE  (DEG) 

MUST  eE  0,15, 30, 45. 60, 75,  OR  90 
IF  90,  THEN  ALL  THE  PREVIOUS  6  LATITUDES 
WILL  BE  COMPUTED 
I MOL  -  MOLECULAR  SPECIE 
I  .  02 

2.  0 

3.  N2 

4.  N 

5.  HE 

6.  AR 

7.  H 

8.  ALL  OF  THE  ABOVE 

IDIFEQ  -  0  DIFFUSIVE  EQUILIBRIUM  NOT  ASSUMED 
NOT  0  DIFFUSIVE  EQUILIBRIUM  ASSUMED 
CARD  3.  COL  1-4,  A  FORMAT 

MMMM  -  SOLAR  ACTIVITY  INDICATOR, 

*GMIN*  OR  "GMAX" 

IF  NEITHER  THEN  BOTH  MIN  ANO  MAX  WILL  BE  DONE 
CARD  4,  COL  1-4,  A  FORMAT 

DODD  -  TIDAL  ZONAL  WAVE  NUMBER  INDICATOR, 

"OIUR"  OR  “SDIU" 

IF  NEITHER  THEN  BOTH  OIUR  ANO  SOIU  WILL  BE  DONE 
CARD  5,  CCL  1-5,  15 

I DUMP 1  -  IF  NOT  ZERO, INTERMEDIATE  DUMP  IS  TRIGGERED 


Figure  1.  Program  TIDEVAR  Inputs  and  Outputs 
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TAPE1  -  CARO  IMAGES 

CAROS  1-25  -  REDUCED  HEIGHTS  AND  GEOMETRIC  HEIGHTS  FOR  SOLAR 
MINIMUM 

XH(N.  1  ) ,2H(N. 1 ) ,n-1 ,25  (3X,F5.2,3X,F5. 1  ) 

XH  *  REDUCED  HEIGHT  (NO.  OF  SCALE  HEIGHTS  FROM  GROUND) 

ZH  .  GEOMETRIC  HEIGHT  (KM) 

CARDS  26-50  “  REDUCED  HEIGHTS  AND  GEOMETRIC  HEIGHTS  FOR  SOLAR 
MAXIMUM,  IN  SAME  FORMAT  AS  FOR  MINIMUM. 

THESE  ARE  FOLLOWED  BY  FOUR  GROUPS  OF  151  CARDS  EACH 
GROUP  1.  DIURNAL  TIDES,  SOLAR  MINIMUM 
GROUP  2.  DIURNAL  TIDES,  SOLAR  MAXIMUM 
GROUP  3.  SEMIDIURNAL  TIDES,  SOLAR  MINIMUM 
GROUP  4.  SEMIDIURNAL  TIDES,  SOLAR  MAXIMUM 
IN  EACH  GROUP  THE  CARDS  ARE  FORMATTED  AS  FOLLOWS 
FIRST  CARD  (CURRENTLY  UNUSED) 

NZONAL , NVODE  ISEN( I ) . 1*1 ,70 
2I2.70A1 

REMAINING  150  CARDS 

(( (U( I ,d.K) ,UH( I , J.K) . I»1 .4) . J«1 ,6) ,K«1 ,25) 

(F6.2.F4. 1 ,F6.2,F4.1 .F6.4.F4.1 .F6.2.F4.1 ) 

U ( I , J , K )  =  AMPLITUDE  OF  ITH  TIDAL  FIELD, 

AT  JTH  LATITUDE  AND  KTH  HEIGHT 
UH ( I , U , K )  *  CORRESPONDING  PHASE  (HOURS) 

TIDAL  FIELDS 

1.  EASTWARD  WIND  VELOCITY  (M/SEC) 

2.  SOUTHWARD  WIND  VELOCITY  (M/SEC) 

3.  UPWARD  WIND  VELOCITY  (M/SEC) 

4.  TEMPERATURE  (K) 

LATITUDES  (DEGREES) 

1  .  0 

2.  15 

3.  30 

4.  45 

5.  60 

6.  75 

HEIGHTS  ARE  GIVEN  IN  THE  XH  AND  ZH  ARRAYS 
-  U70  BACKGROUND  DENSITIES 

BINARY  RECOROS.  THE  FIRST  FOR  SOLAR  MINIMUM, 

SECOND  FOR  SOLAR  MAXIMUM 
( HBACK ( J , I ) , (BACKGR(d,K. I ) ,K-1 , 6 ) , RHO( d , I ) , 

J»1 .23) 

HBACGR  -  HEIGHT ( KM) 

BACKGR  -  BACKGROUND  NUMBER  DENSITIES  (LOG  BASE  10, 

MXS) 

1  .  N2 

2.  02 

3.  0 

4.  A 

5.  HE 
6  •  H 

RHO  -  TOTAL  MASS  DENSITY  (MKS) 

OUTPUTS 

FOR  EACH  SOLAR  ACTIVITY  LEVEL  AND  ZONAL  WAVE  NUMBER  REQUESTED. 

THE  TOTAL  DF.NSITY  TIDE  IS  PRINTED  OUT  AND  WRITTEN  TO  TAPE  4  IN 
THE  BINARY  FORMAT 

L , (HH( KK) , ( RAMP( JJ  ,  KK ) , RPHASE ( JJ , KK) , JJ«1 , 6 ) ,KK. 1 , L) 

WHERE 

L  *  NUMBER  OF  HEIGHTS 

HH  «  GEOMETRIC  HEIGHTS  (KM) 

RAMP ( J J , KK  )  «  AMPLITUDE  AT  HEIGHT  HH( KK ) , LAT I TUDE  (JJ-1)-15  DEGREES 
RPHASE(dd.KK)  *  PHASE  IN  HOURS,  AT  THAT  HEIGHT  AND  LATITUDE 

IF  I DUMP 1  ■  1,  DETAILED  COMPOSITIONAL  TIDES  ARE  ALSO 
PRINTED  OUT 


TAPE3 

TWO 

THE 
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makes  a  negligible  contribution  to  the  total  mass  density,  as  does  H  in 
the  altitude  ranges  considered  (h^x  =  283km  for  solar  minimum,  420km 
for  solar  maximum).  The  computations  can  be  carried  out  for  6 
latitudes  0,  15,  30,  45,  60,  and  75  deg.,  for  both  minimum  and  maximum 
solar  activity  (corresponding,  approximately,  to  exospheric 
temperatures  800K  and  1400K,  respectively),  and  for  both  diurnal  and 
semidiurnal  tides.  Results  are  combined  with  appropriate  time  averaged 
background  distributions  to  yield  total  mass  density  tides  which  are 
written  to  disk  for  later  use. 

The  diffusion  coefficients  are  computed  as^ 

t  b. 

ai  Tn  1 

Di  =  (2705) 

where 

N0  =  total  time  averaged  number  density 
T0  =  time  averaged  temperature 

and  a-j  and  b-j  are  constants  given  in  Table  1,  along  with  the  thermal 
diffusion  coefficients  <*•$. 

TABLE  1.  Thermal  and  Molecular  Diffusion  Parameters 


Gas 

“i 

ai  (m-i.  s-1) 

bi 

02 

0 

4.863x1020 

.75 

0 

0 

6.986x1020 

.75 

n2 

0 

6.986x1020 

.75 

N 

0 

7.5x1020 

.75 

He 

-0.38 

1.7x1021 

.691 

A  r 

0 

4.487x1020 

.870 

H 

-0.25 

3.305X1021 

.5 
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The  background  atmosphere,  defining  N0  and  T0  ,  is  a  latitude 
independent  model  defined  by  Hong  and  Lindzen'^. 

Time  averaged  vertical  flux,  important  only  for  H  and  chemical  loss, 
have  been  neglected,  although  they  can  be  easily  incorporated  later,  if 
warranted,  for  detailed  composition  studies.  In  the  altitude  range  of 
interest,  H  is  a  minor  specie;  and  the  only  important  chemical  loss 
process,  0+  to  O2  charge  transfer, ^  affects  only  the  diurnal 
variation  of  O2.  Thus  neither  of  these  should  have  an  important  effect 
on  total  mass  density  semidiurnal  tides. 

The  tidal  perturbations  T',  V,  and  W,  are  those  generated  by  Forbes  and 
Garrett16  for  the  diurnal  tide  and  Garrett  and  Forbes11  for  the 
semidiurnal  tide. 

The  numerical  integration  is  performed  with  respect  to  the  reduced 
height  variable, 
y  —  fh  dh 

A  "  Jo  T 

where  H  is  the  mean  molecular  scale  height  according  to  the  background 
model.  With  this  change  of  variable  the  various  derivatives  in  the 
governing  equation  for  Rj  transform  as: 


JL  .  _i_  ax 
3h  ~  3x  3h 


1  J_ 
H  3X 


.  i  (1  _3_  \  «_1  3H  J_  +1  3 2 
3x  H  3X  '  H2  3X  3x  H  Fit2 


In  terms  of  the  new  variable  x,  it  is  feasible  to  use  a  uniform 
numerical  integration  grid. 


1 


The  upper  boundary  conditions  at  h  s  283km  and  421km  for  minimum  and 
r aximum  solar  activity,  respectively,  are  defined  by  diffusive 
equilibrium: 


At  the  lower  boundary,  100km,  an  attempt  was  made  to  calibrate  the 
value  of  using  satellite-borne  mass  spectrometer  data  at  140 
-160km^»18.  This  feature  is  implemented  aa  "Option  2"  in  the  code,  in 
contrast  to  the  normal  "Option  1”  defined  by  preset  values  of  at 
100km,  and  the  diffusive  equilibrium  condition,  given  above,  at  the 
upper  boundary.  This  effort  proved  fruitles,  however,  as  the  solutions 
at  140- 160km  were  quite  insensitive  to  the  assumed  boundary  condition. 
This  ia  evidently  due  to  the  smallness,  at  100km,  of  molecular 
diffusion,  ttfiich  is  the  only  process,  in  the  model  being  used,  by  which 
molecules  at  100km  can  move  vertically  upward  to  affect  distributions 
at  higher  altitudes.  It  is  possible  that  the  inclusion  of  eddy 
diffusion  in  the  model  could  change  this  situtation.  However,  based  on 
Reference  14,  eddy  dif fusivities  exceed  molecular  diffusivities  by  at 
most  a  factor  of  2  above  100  km  for  the  most  important  constituents. 
Therefore  the  attempt  to  calibrate  the  lower  boundary  conditidn  was 
dropped,  and  the  fixed  condition 

Ri  =  0 

was  adopted. 

The  numerical  integration  method  employed  is  that  described  by  Lindzen 
and  Kuo^  (for  one  uncoupled  differential  equation)  and  Richtmeyer^O . 

The  Jacchia  197021  density  model  has  been  used  to  specify  the 
background  composition,  since  this  has  been  recommended  as  the  best 
model  for  satellite  use22 ,  The  more  recent  Jacchia-Bass  and  MSIS 
models  are  more  expensive  computationally  than  the  Jacchia  1970  model 
while  yielding  no  significant  improvements  in  accuracy^.  The  Jacchia 
1971 24  model,  similar  to  the  1970  model  in  expense  and  total  density 
accuracy,  represents  composition  less  accurately  than  the  Jacchia  1970 
model2* . 
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The  resulting  semidiurnal  tide  was  fit  to  the  following 
represent  at  ions : 

Region  1:  120km  -  145km 

Alnd  =  X  (  h,  e)  cos2at  +  Y(  h,  9)  sin2at 


where 


=  total  mass  density 
=  altitude  (km) 

=  latitude 

=  earth's  rotation  rate 
=  local  t ime 


2  2 


X  «  co«£  6  -lQ  ( A(1  ,j)  +  B(i  ,j)  (TINF  -  800)} 


• (h-145)1  sin2j  e 


2  2 


Y  =  cos  6  i|Q  jIq  (C(i,j)  +  D(1  ,j) (TINF  -  800 J) 


•(h-145)1  sin2je 


TINF  s  global  mean,  geomagnet ical ly  quiet,  exospheric 
temperature  (K) 

A(i»j),  B(i,j),  C(i,j),  and  0(i,j)  s  adjustable  parameters 
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Region  2:  155km  -  h^.  h^,  =  283km  -«-  0.23  (T INF -800) 


Same  as  Region  1  except: 


3  2 


X  =  cos  e  ^  J.|0  {E (i  ,j)  +  F(1,j)  (TINF-800)} 

•  (h-155 ) 1  sin^e 


3  2 

Y  =  COS2  6^  {6(i, j)  +  H(i,j)  (TINF  -  800)} 

^  (h-155)1  sin2je 
Region  3:  145km  -  155km 

Cubic  polynomial  which  matches  Region  1  function  and  slope  at  145km  and 
Region  2  function  and  slope  at  155km. 


Region  4:  h  greater  than  ly, 


Alnd  =  FAjlnd  +  (l-FjAglnd 


where 


F  *  M0no  exp  {-b0(h-hm))/KT 

Hj  «  M0no  exp  {-b0(h-hm)}+  exp  t -bHfi( h-hj } 


Mq  -  16;  MHe  =  4 


log  n0(/cm3)  =  8.72-0.00063  (TINF-800) 
log  nHe(/cm3)  -  6.84-0.0003  (TINF-800) 
b1  -  M1g/(R-TINF) 
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g  =  9.8m/sec 

R  *  8314.32  J(Kg-Mo1 )" Vk 

3Alnd(h  ) 

ijlnd  ■  AlndOg  +  (h-hj  — m  • 

A2lnd  =  lndQ  (h,  TINF  +  aTINF)  -  lndQ(h,  TINF) 

d0  =  daily  average  static  diffusion  density,  frorrTliabTes 
ATINF  =  cos2  8  (c  +  b  sin20) 

c  *  cos2at  <6.4  +  0.086  (TINF-800)}  +  sin2at  (61.15  +  0.0175  (TINF800)} 
b  =  cos2at  (97.3-0.221  (TINF-800)}  +  sin2at  (-156.8  +  0.23  (TINF-800)} 


Table  2  contains  the  coefficients  for  regions  1  and  2.  It  was  found 
possible  to  fit  these  two  regions  separately  to  low  order  polynomials. 
Fitting  the  combined  region  above  120km  proved  difficult  even  with  10th 
order  polynomials.  Furthermore  a  piecewise  low  order  fit,  as  presented 
here,  is  preferable  to  a  single  high  order  fit  for  practical  use.  The 
latitudinal  dependence  hopefully  simulates  the  first  2  or  3  semidiurnal 
tidal  modes  and  reflects  a  north-south  symmetry.  Inclusion  of 
seasonal-latitudinal  effects  would  require  terms  with  odd  powers  in 
sine  and  functional  dependence  on  the  day  of  the  year. 

The  temperature  dependence  reflects  the  assumption  of  linear  behavior 
between  the  two  limiting  cases  of  minimum  solar  activity  (TINF  =  800K) 
and  maximum  solar  activity  (TINF=1400K) .  The  fits  obtained  in  Regions 
1  and  2  were  quite  good  for  both  minimum  and  maximum  solar  activity 
with  residuals  in  Alnd  generally  less  than  0.03. 

High-latitude  data,  above  45  degrees,  was  excluded  from  the  fits.  An 
IBM  polynomial  regression  program^  was  adapted  for  use  on  this  effort. 
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Table  2.  Coefficients  for  Total  Density  Semidiurnal  Tide  in  Regions  1  and  2 
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The  Region  4  formulation  is  an  attempt  to  extend  the  computations  above 
the  upper  limit,  hm,  at  which  the  wind  and  temperature  tides  have 
reached  their  asymptotic  values.  Above  this  point  the  density  tide 
might  be  expected  to  grow  linearly  in  height,  controlled  by  diffusive 
eguilibrium.  However,  this  condition  does  not  occur  because  there  is 
a  change  in  composition  from  predominantly  atomic  oxygen  at  hm  to 
helium,  and  at  still  higher  heights,  to  hydrogen.  In  the  exospheric 
limit  the  simple  solution,  in  which  the  density  is  controlled  by  the 
temperature,  has  been  adopted.  The  semidiurnal  variation  in  the 
exospheric  temperature  was  obtained  by  2-point  fits  to  the 
Garrett-Forbes^  results  at  0  and  45  degrees. 
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By  comparison  the  original  Jacchia  1970  formulation  is: 

Tj  *  T.  (1  +  Rsin"x)  u  *  r£2 s.Vstffy  cosn  ta  , 

1  +  Rsinmx  T 

where 

n  =  3.0 

ta  =  function  of  local  time 

Thus  the  only  changes  are  in  the  local  time  dependent  term,  the  second 
term  inside  the  second  pair  of  parentheses.  The  local  time  term  in  the 
original  model  is  replaced  by  a  purely  diurnal  term  plus  a  constant 
which  equals  the  average  of  the  original  local  time  term.  The  diurnal 
amplitude  factor  b  is  selected  to  simulate  recent  data1^  and  theory.16 


In  particular,  measured  diurnal  density  variations  at  low  solar 
activity  indicate  that  the  Jacchia  1971  model  diurnal  tide,  which  is 
nearly  the  same  as  in  the  Jacchia  1970  model,  is  too  high.  The 
temperature  amplitude  predicted  by  the  Jacchia  1970  model  is  also 
higher  at  low  solar  activity  than  that  predicted  theoretically,16  while 
there  is  closer  agreement  at  high  solar  activity.  The  phase  depends 
on  the  height  at  which  the  density  is  being  computed.  This  phase 
variation  also  relects  both  data  and  theory  at  low  altitude. 


Eventually,  as  resources  and  time  permit,  an  improved  representation  of 
the  diurnal  density  tide  should  be  constructed,  as  was  done  above  for 
the  semidiurnal  tide,  basing  the  model  directly  on  density,  rather  than 
temperature,  variations.  The  diurnal  representation  presented  here 
should  nevertheless  be  an  improvement  over  that  given  in  the  Jacchia 
1970  model,  and  would  be  most  inaccurate  at  low  altitudes,  where  the 
semidiurnal  tide  dominates. 


1.1.2  3acchia-Bass  Tidal  Model 


In  this  subsection  we  discuss  the  development  of  an  improved  tidal 
model  for  the  3acchia-Bass  density  model^.  The  Jacchia-Bass  (3B)  model 
is  an  adaptation  of  the  3acchia  1977  model, 26  incorporat inq  an 
analytically  integrable  temperature  profile  for  the  static  diffusion 
equation  above  125  km.  This  feature  significantly  reduces  the  computer 
memory  requirements  for  the  solutions  to  the  static  diffusion  equation 
(which  in  the  original  model  must  be  tabulated  numerically  for  each  of 
7  components)  and  thus  significantly  broadens  the  range  of  computers  on 
which  the  model  can  be  run. 

The  tidal  representation  in  the  3acchia  1977  model,  and  hence  also  in 
the  3B  model,  is  significantly  deficient  in  the  semidiurnal  tide  at  low 
altitudes,  as  indicated  in  the  previous  subsection.  In  an  attempt  to 
capture  some  detail  in  the  diurnal  tides  not  accommodated  in  previous 
models,  a  height -dependent  exospheric  temperature  diurnal  phase  is  used 
in  the  computation  of  the  number  density  of  each  constituent. 

The  model  presented  here  includes  a  model  of  the  temperature  diurnal 
and  semidiurnal  tides  above  125  km  which  retains  the  analytic 
integrabil ity  of  the  3B  model.  Analytic  expressions  for  the  component 
number  density  tides  are  then  immediately  obtained,  with  the  exception 
of  boundary  conditions  and  corrections  for  the  effects  of  winds.  Winds 
should  be  important  above  125  km  only  for  atomic  oxygen  and  minor 
constituents.  In  the  upper  thermospheric  region  where  static  diffusion 
is  valid  (wind  effects  negligible),  boundary  conditions  can  be  obtained 
from  the  solutions  to  the  equations  governinq  the  composition  tidal 
perturbation  discussed  in  the  previous  subsection.  Departure  from 
static  diffusion  at  lower  altitudes  can  then  be  obtained  by  comparison 
of  the  exact  and  static  diffusion  solutions. 

Below  125  km,  where  departures  from  static  diffusion  are  large,  it  is 
recommended  that  direct  fit  to  the  exact  Forbes  composition  tides  be 
employed. 
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The  temperature  model  that  has  been  developed  above  125  km  is  as 
follows: 


(A„  cosmut  +  B  sinmut) 
m  m 


where 


T  =  temperature 

a)  =  rotation  rate  of  earth 

t  =  local  time 


m 
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_  p_cmx 
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q  m 
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£ 

1=0 


mi 


x  =  geopotential  height  above  125  km 

am’  cm*  dmf  ’  bm’  fmj  are  funct*ons  of  solar  activity  and  latitude. 


The  equinox  temperature  tides  computed  by  Forbes  and  Garrett  (see 
previous  subsection)  have  been  used  to  determine  the  parameters  am,  bm, 
cm>  ^mi*  *mi  the  latitudes  0,  15,  30,  45,  60,  and  75  degrees  and  for 
minimum  and  maximum  solar  activity  (F10.7  =  75  and  260,  respectively, 
where  F10.7  is  the  10.7  cm  solar  flux  in  10"22  watts/m2/Hz) . 

Least  squares  multiple  regression  program  STEPR25  was  adapted  for  this 
purpose.  In  this  program  Laguerre  polynomials2^  were  used  in  place  of 
powers  of  the  geopotential  height.  This  was  done  in  the  interest  of 
any  efficiency  that  might  be  gained  in  the  fitting  process  from  the  use 
of  orthogonal  functions.  A  post-processing  program  LAGTP0W  converts 
the  expansions  in  Laguerre  polymomials  to  power  series.  An  iterative 
procedure  was  added  to  find  the  optimum  parameter  cm,  for  each  latitude 
and  solar  activity  level.  For  fixed  value  of  cm  a  regression  is 
performed  to  determine  the  other  parameters.  Then  a  one  dimensional 
search  over  cm  is  made  to  find  the  value  for  which  the  sum  of  the 
uquared  residuals  minimizes.  Satisfactory  fits  were  obtained  by 
truncating  the  power  series  in  the  geopotential  height  x  at  n=1. 


The  parameters  am,  bm,  dm^,  and  fm^  are  next  fit  as  functions  of 
latitude  to  expansions  in  the  associated  Legendre  functions  P  .  Only 
even  n  is  used  since  symmetry  about  the  equator  is  assumed.  The 
nonlinear  parameters  cm  are  fit  to  expansions  in  Legendre  polynomials, 
again  including  only  the  even  n  functions.  Program  LATFIT  was 
constructed  from  program  STEPR  to  perform  these  latitudinal  functional 
determinations. 

Given  the  temperature  profiles  above  it  is  possible  to  obtain 
immediately  the  expressions  for  the  tidal  perturbations  in  the  number 
densities  under  the  static  diffusion  approximation: 

MnNK  =  GK  +  (l+aK)Aln  |  -MKAF(h) 

where 

=  number  density  of  constituent  K 
h  =  altitude 
Gk  *  Ain  Nk(h0) 
hQ  =  reference  altitude 

«K  =  thermal  diffusion  coefficient  of  constituent  K 
T  =  temperature 

M.  =  molecular  weight  of  constituent  K 
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Updates  to  Programs  CADNIP  and  BADMEP 

Programs  CADNIP  and  BADMEP^  have  long  been  used  by  AFGL  scientists  in 
atmospheric  density  modeling.  Program  CADNIP  is  used  for  density 
determination  by  finding,  in  a  least  squares  sense,  the  scale  factor 
which,  when  multiplied  by  a  chosen  density  model,  yields  by  numerical 
integration  the  resulting  satellite  ephemeris  which  best  fits  available 
tracking  data.  BADMEP  evaluates  models  by  comparing  ephemerides 
generated  with  them  to  available  tracking  data.  This  section  will 
discuss  the  following  updates: 

o  Improvement  of  routines  for  converting  between  mean  elements 
and  postion/velocity  vectors 

o  Addition  of  options  to  input  and  output  position/velocity 
directly 

o  Temporary  modifications  to  test  variants  of  an  existing 
density  model 

Discussion  of  these  will  be  followed  by  addenda  to  the  CADNIP  and 
BADMEP  User's  Guides  (Appendices  A  and  B  of  Reference  2)  bringing  the 
user  fully  up  to  date  on  usage  of  the  two  programs. 

1.2.1  Mean  Elements  -  Position/Velocity  Transformation 

Until  recently  CADNIP  and  BADMEP  have  required  input  of  mean  Kepler ian 
elements  at  some  specified  time  to  initiate  processing.  In  CADNIP  they 
are  used  to  develop  an  initial  estimate  of  the  ephemeris  for  the 
iterative  least  squares  analysis.  In  BADMEP  they  are  used  to  initiate 
the  ephemeris  generation  for  evaluation  of  the  density  models.  In 
either  case  the  elements  must  first  be  transformed  to  a 
posit  ion/velocity  vector  since  the  numerical  integration  is  done  in 
cartesian  coordinates.  In  addition  CADNIP  also  outputs  the  mean 
elements  at  epoch  for  the  best  fit  orbit,  which  requires  the  reverse 
transformation:  from  position/velocity  to  mean  elements. 


The  transformation  from  mean  elements  to  position/velocity  requires  the 
following  two  steps: 

1)  Add  to  the  elements  the  short  periodic  corrections  due  to  the 
second  harmonic  32  the  geopotential 

2)  Transform  the  resulting  osculating  elements  to 
pos it ion/veloc it  y 

For  the  inverse  transformation  the  steps  are  done  in  reverse:  the 
osculating  elements  are  formed  from  the  posit  ion/ velocity  vector,  and 
from  these  the  short  periodic  corrections  are  subtracted.  Note  that 
this  definition  of  the  transformation  implies  that  the  mean  elements 
are  "mean"  only  with  respect  to  the  short-periodic  perturbations.  This 
may  differ  from  the  "mean"  elements  supplied  by  some  outside  agencies 
such  as  ADCOM,  whose  "mean"  elements  also  exclude  long-periodic 
perturbations  due  to  the  third  harmonic  33.  Such  elements  should  be 
used  with  care.  However  in  the  most  common  use  of  such  element  sets  in 
the  CADNIP/BADMEP  system,  as  initial  estimates  to  be  refined  by  CADNIP, 
extreme  accuracy  is  not  a  neccessity. 

A  large  source  of  error  arises  in  the  computation  of  the  3hort -per iodic 
corrections  to  the  elements  for  nearly  circular  orbits,  due  to  a 
singularity  caused  by  loss  of  perigee  definition.  This  has  resulted  in 
such  anomalies  as  "negative"  eccentricities  appearing  in  the  printouts. 
Aksnes2®  has  developed  a  formulation  in  which  the  perturbations  are 
computed  for  a  set  of  intermediate  coordinates,  the  Hill  variables: 

r  =  radius  vector 

r  =  time  derivative  of  r 

G  =  ka(1-e2) 

H  =  G  cos  i 

u  =  g  +  f 

h  =  right  ascension  of  ascending  node 


k  =  earth's  gravitational  constant 
a  =  semimajor  axis 
e  =  eccentricity 
i  =  inclination 
g  =  argument  of  perigee 
f  =  true  anomaly 

The  Hill  variables  remain  well  defined  for  circular  orbits  and  thus  the 
perturbations  remain  non-singular. 

To  enable  the  treatment  of  nearly  circular  orbits,  the  subroutines  in 
CADNIP  and  BADMEP  performing  this  transformation  have  been  replaced 
with  the  following  routines. 
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OSCTMN: 

MNTOSC: 


PVTHIL: 

HILTMN: 

ECCF: 


MNTHIL: 

SHP: 

HILTPV: 


main  subroutine  to  convert  position  and  velocity  to  mean 
elements 

main  subroutine  to  convert  mean  elements  to  position  and 
velocity 

converts  position  and  velocity  to  Hill  variables 
converts  Hill  variables  to  mean  elements 
computes  true  anomaly  f,  ecosf,  esinf,  and  mean  anomaly 
M,  given  r,  r  and  G 

converts  mean  elements  to  Hill  variables 

computes  short-periodic  corrections  to  Hill  variables 

converts  Hill  variables  to  position  and  velocity 


In  the  computation  of  the  Hill  coordinates  from  posit  ion/velocity  we 
have 

r  =  /  x2  +  yZ  + 
r  -  (xx  +  yy  +  zz)/r 


Gx  -  /K  (yz  -  zy) 


Gy  s  /IT  (zx  -  xz) 


g7  =  /TT  (xy  -  yx) 


•  .*  o  ; 


sin  u  =  G  £in  («|>-h)  cose/H 
cos  u  =  cos  (4>-h)  cose 

where 

0  =  latitude 

Note  that  in  calculating  u,  the  cose  factor  cancels  out  (sinu/cosu), 
and  thus  need  not  be  calculated.  Exceptions  are  easily  handled  for 
equatorial  and  polar  orbits.  In  the  case  of  equatorial  orbits,  h  is 
arbitrary,  and  therefore  set  to  0,  with 

u  =  ♦ 


For  polar  orbits  we  have 

sin  u  =  z/r 

cos  u  =  SGN(z)  /1-sin^u 

where  SGN(X)  is  the  algebraic  9ign  of  X. 


To  compute  the  mean  elements  from  the  Hill  variables,  we  compute 
e  sinf  =  rG/k 

e  cosf  *  G2/ (kr)-l 

r  _  1-e2 
a  1  +  e  cosf 

e  sinE  = 

e  cos  E  =  1  -- 

d 

* 

M  =  £  -  esin  E 
g  =  u  -  f 
i  =  cos-1  (H/G) 
a  =  { G2/ ( i -e2)}  /k 

To  compute  the  Hill  variables  from  the  mean  elements,  first  the 
eccentric  anomaly  is  E  obtained  by  solving  Kepler's  equation 

M  =  E  -  esin  E 


resinf 

a/l-e^ 


Then  we  have 

r  =  a(l  -  e  cosE) 

cosf  =  (cosE-e)/(l-e  cosE) 

sinf  =  (  *{l-e2  sin  E)/(l-e  cosE) 

u  =  f  +  g 

G  =  /ka(l-e2) 

H  =  G  cosi 

r  *  k  esinf/G 
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To  compute  the  position  and  velocity  from  the  Hill  variables  requires 


xm  ’  -6  s1nh 

ym-  |  cost, 


x  =  cosh 
n 

yn  *  sinh 

2n'° 
s  *  G/r 

x  *  r  Umsin  u  +  xn  cos  u) 

y  =  r  (y  sin  u  +  y„  cos  u) 

m  n 

z  =  r  (zmsin  u  +  z„  cos  u) 
m  n 

x  *  rx/r  +  s(xm  cos  u  -  xp  sin  u) 

y  =  ry/r  +  s(ym  cos  u  -  yn  sin  u) 

z  =  rz/r  +  s(zm  cos  u  -  z„  sin  u) 

m  n 

The  computation  of  the  short -per iodic  corrections  to  the  Hill  variables 
is  done  in  accordance  with  Reference  28,  including  only  terms  due  to 
the  second  harmonic  32. 
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A  significant  portion  of  the  internal  computations  of  CADNIP  and  BADMEP 
is  performed  in  a  system  of  canonical  units  of  length  and  time, 
specified  by: 

o  unit  of  length  =  earth  radius 

o  unit  of  time  =  time  required  for  a  satellite  in  an  earth 

radius  circular  orbit  to  travel  1  radian. 

In  these  units  the  earth's  gravitational  constant  k  is  unity  and  thus 
will  not  appear  in  equations  expressed  in  these  units.  Thus  one  will 
notice  the  absence  of  k  on  examining  much  of  the  CADNIP/BADMEP  code, 
including  that  just  discussed. 

1.2.2  Input  and  Output  of  Position/Velocity 


Options  have  been  added  to  CADNIP  and  BADMEP  to  input  position  and 
velocity  instead  of  elements  and  to  output  either  or  both.  In  the  case 
of  CADNIP  the  user  indicates  by  an  input  flag  whether  the  starting 
estimate  of  the  satellite  ephemeris  is  a  position/velocity  vector  or  an 
element  set.  At  the  end  of  each  iteration  in  the  differential 
correction  procedure,  if  such  printout  is  desired,  and  at  the  end  of 
the  final  iteration,  CADNIP  will  print  the  posit  ion/velocity,  and/or 
the  element  set,  at  the  epoch  of  the  fit  span,  in  accordance  with  a 
second  input  flag.  The  units  in  the  printouts  are  km  and  km/sec,  and 
deg.  Punched  cards,  in  NAMELIST  format  appropriate  for  BADMEP,  are 
produced,  if  requested,  for  the  final  values  of  both  elements  and 
position/velocity  at  the  epoch  of  the  fit  span.  For  an  element  set  the 
NAMELIST  name  is  NEWIN,  as  previously,  but  for  a  position/velocity 
vector  the  NAMELIST  name  is  NEWINX.  The  punched  position/velocity 
vector  is  in  canonical  units  (see  end  of  Section  1.2.1). 

BADMEP  will  read  its  run  parameters  from  the  NAMELIST  group  NEWIN  or 
the  group  NEWINX  according  to  an  input  flag  specified  on  the  preceding 
device  control  card. 
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The  two  group 8  differ  only  in  that  NEWIN  includes  an  element  set,  while 
NEWINX  includes  a  position/velocity  vector.  The  cards  punched  by 
CADNIP  include  the  time  of  the  elements  or  vector  and  the  model  density 
factor,  as  well  as  the  element  set  or  vector.  Other  inputs  must  be 
inserted  by  the  user. 

Details  of  the  operation  of  these  options  and  other  features  are 
described  in  Section  1.2.4. 


1.2.3 


Thermospheric  Response  to  High  Solar  Activity 


AFGL  scientists  conducted  orbital  drag  studies  over  a  period  of 
extremely  high  solar  activity,  4-16  November,  1979,  with  the  objective 
of  determining  a  suitable  modification  to  existing  empirical  models 
for  such  high  activity.  The  Jacchia  197021  model  was  chosen  for  this 
study.  Several  modifications  of  the  formula  for  the  global  nighttime 
minimum  exospheric  temperature  at  zero  geomagnetic  activity  were 
developed  by  AFGL  researchers.  These  modifications  were  implemented  in 
software  and  made  available  to  AFGL  scientists  for  use  in  these 
studies.  The  modifications  were  confined  to  subroutine  APFTMP  which 
tabulates  solar  flux  and  geomagnetic  activity  dependent  quantities  as 
functions  of  time.  Details  of  these  studies  are  given  elsewhere .29 

1 .2 .4  Addenda  to  CADNIP  and  BADMEP  User's  Manuals 

The  following  are  addenda  to  the  CADNIP  and  BADMEP  User's  Manuals  which 
are  provided  in  Appendices  A  and  B  of  Reference  2. 

1.2 .4.1  Density  Model  Preparation 


The  Jacchia  1977  model  requires  a  special  mass  storage  file  "JSDM" 
whose  format  is  given  in  Table  3. 

1.2 .4 .2  Solar  and  Geophysical  Data 


The  format  of  the  solar  and  geophysical  file  is  still  as  described  in 
Section  A.4,  Reference  2.  However  the  data  must  begin  at  least  2  days 


prior  to  the  start  of  the  period  to  be  processed  (3  days  if  the  Jacchia 
1977  model  is  used),  and  there  must  be  a  minimum  of  308  days  of  data 
(462  if  the  Jacchia  1977  model  is  used). 


1.2. 4. 3  CADNIP  Starting  Elements  Cards 

As  indicated  in  Section  1.2.2,  either  starting  elements  or  starting 
position  and  velocity  may  be  supplied.  The  user  must  indicate  which  by 
an  input  option  flag  in  columns  73-75  of  the  first  card  of  this  2-  card 
set.  The  value  of  this  flag  must  be  1  for  position  and  velocity,  any 
other  value  for  elements.  An  output  option  flag,  in  columns  76-78  of 
the  first  card,  indicates  whether  elements,  position  and  velocity,  or 
both,  at  the  epoch  of  the  fit  span,  are  to  be  printed  out  in  the 
results  of  the  differential  correction  procedure. 

Flag  Value  Printout  Includes 

0  elements  only 

1  position  and  velocity  only 

any  other  value  both 

The  remaining  parameters  of  the  first  card  are  given  as  on  page  25, 

Reference  2,  exept  substitute  "position  and  velocity"  for  "elements"  if 
appropriate. 


If  position  and  velocity  are  to  be  input,  they  must  appear  on  the 
second  card  as  follows: 


Columns 


Description 


1-10 

11-20 

21-30 

31-40 

41-50 

51-60 


x  position 
y  position 
z  position 
x  velocity 
y  velocity 
z  velocity 


(km) 

(km) 

(km) 

(km/sec) 
( km/sec ) 
( km/sec ) 


In  this  case  the  card  is  read  with  the  format  (3F10.4, 


3F10.5). 


Table  3.  Format  of  Mass  Storage  File  "JSDM" 
for  CAONIP  and  BADMEP 


Record  Number  1 


Word 


Symbol 


Description 


1 

LABEL 

2-9 

TITLE 

10 

TMINMOD 

11 

TMAXMOD 

12 

NTSTEP 

13 

DTMOD 

14 

NREG  (<7) 

15 

NSUM 

16-23 

ALZSTEP^ 

24-30 

NZSTEPi 

31-37 

ALDZi 

38 

GNMOD 

39 

RNMOD 

40 

RSTAMOD 

41 

AVOGMOD 

42 

NSPMOD  (6) 

43-48 

AMWM0Di 

49-55 

M0DRLENi 

[“(NSPM0D+1) 

♦(NTSTEP+1) 

*(NZSTEPi+l)] 

The  species  will  be  given  in  the  order: 


BCD  label 
80  Character  BCD 
Minimum  exospheric 
temp,  in  table 
Maximum  exospheric 
temp. 

Number  of  Temperature 
steps 

Temperature  step  size 
Number  of  height  regions 
(=number  of  subsequent 
records ) 

Unused 

Natural  logs  of  boundaries 
(km)  of  height  regions. 
Number  of  height  steps  in 
each  region 

Step  size  in  natural  log 
of  height  in  each  region. 
Gravitational  acceleration 
at  Earth  surface  (km/sec2) 
Earth  radius  (km) 

Universal  gas  constant 
Avogadro’s  number 
Number  of  species 
Molecular  masses  of 
species 

Lengths  of  subsequent 
records 


02,  0,  N2,  HE  AR,  H 


Pj| 
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Table  3  (Cont'd. ) 


Record  Number 

Word 

1 

2  thru 
NSPMOD+1 

NSPMOD+2 

thru 

(NTSTEP+1)* 

(NSPMOD+1) 

(NTSTEP+1 )*(  'NSPMOD+1 )+l 
thru 

MODRLEN  (I REG) 


IREG+1  1<IREG<NREG) 


Description 

Local  temperature  at  height  exp  [ALZSTEP(IREG)] 
and  exospheric  temperature  TMINMOD 
Logs  base  10  of  number  densities  (m_3)  for  molecular 
species  at  this  height  and  exospheric  temperature. 

Repeat  of  words  1  thru  NSPMOD+1  for  this  height  and 
equally  spaced  exospheric  temperatures  TMINNOD  + 
DTMOD  thru  TMAXMOD 


Repeat  of  words  1  thru  (NTSTEP+1 )*(NSPM0D+1)  for 
remaining  heights 
exp  [ALZSTEP  (I REG)] 

+  I  *ALDZ  (IREG)] 

I  =  1  thru  NZSTEP  (IREG) 


1.2. 4. 4  CADNIP  Change  Cards 


The  number  of  the  parameter  to  be  changed  (field  2,  page  26, 
Reference  2)  should  appear  in  columns  16-19.  The  increase  in  the  field 
size  to  4  columns  accommodates  the  increase  in  the  number  of  parameters 
(Section  1.2. 4. 6). 

1.2. 4. 5  CAONIP  Run  Card 

The  currently  available  density  models  and  the  indicators  (field  12) 
are: 

1  -  Jacchia  1964 

2  -  1966  Supplements 

3  -  Jacchia  1971 

4  -  USSR  Cosmos38 

5  -  Jacchia-Walker -Bruce 

6  -  Jacchia  1977^6 

7  -  Lockheed/NASA 

8  -  MSIS  778*9 

9  -  1962  U.S.  Standard 

10  -  MSIS  7831 

11  -  DENSEL 

12  -  Jacchia  1970 

13  -  Jacchia  1973 

14  -  Forbes-Garrett-Gillette  Model  B32 

All  the  models  except  4,  6,  8,  9,  10,  11,  and  14  reguire  density  tables 
as  described  in  Sect  ion  A.  3,  Reference  2.  Model  6  reguires  the  mass 
storage  file  JSDM,  as  described  in  Section  1.2. 4.1  of  this  report.  All 
the  models  except  9  and  11  require  solar  and  geomagnetic  activity 
indices  as  described  in  Section  A. 4,  Reference  2.  Those  models  without 
footnotes  in  the  above  list  are  described  in  Reference  2,  Appendix  C. 


1.2. 4. 6  Description  of  CADNIP  Parameter  Table 


The  number  of  parameters  has  been  increased  to  1165  to  accommodate  a 
gravity  model  with  maximum  m,n  =  32  (previous  maximum  was  20).  The 
gravity  model  coefficients  are  arranged  as: 

Cnm:  Parameter  number  80  +  n  (n  +  1)/2  +  m  -  2 

Sp,,,,:  Parameter  number  638  +n  (n  -  1  )/2  +  m  -  1 

Parameter  number  13,  if  not  zero,  suppresses  printout  of  the  gravity 
model  coefficients,  parameters  81  -  1165.  The  default  value  of  this 
parameter  is  0. 

If  parameter  number  49  is  greater  than  zero,  only  the  preliminary 
adjustment  procedure  is  run,  the  differential  correction  procedure  is 
skipped,  and  a  disk  file,  TAPE10,  is  written  of  all  the  tracking 
observations  not  rejected  during  the  preliminary  adjustment  procedure. 
This  file  is  useful  for  subsequent  runs  of  BADMEP,  which  does  not 
filter  out  bad  observations.  The  default  value  of  this  parameter  is  0. 

If  parameter  number  53  is  greater  than  zero,  the  final  values  of  the 
elements  and  position/velocity  at  the  fit  span  epoch  will  be  punched 
out  in  the  NAMELIST  format,  group  name  NEWIN  and  NEWINX,  respectively, 
compatible  with  input  required  by  BADMEP.  The  default  value  of  this 
parameter  is  1 . 

1.2. 4. 7  Output  from  CADNIP 

All  specification  cards  except  the  change  cards  are  listed  in  the 
printout.  This  is  followed  by  a  printout  of  the  non-zero  values  of  the 
parameter  table,  except  that  parameters  81-1165  are  omitted  if 
parameter  number  13  is  not  zero.  This  is  followed  by  the  results  of 
the  preliminary  adjustment  procedure,  and  the  differential  correction 
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procedure.  The  latter  includes  position  and  velocity  if  requested  in 
accordance  with  Section  1.2. 4. 3.  The  elements  and  postion/velocity  are 
punched  in  NAMELIST  format,  as  described  in  Section  1.2.2,  if  parameter 
53  is  greater  than  zero. 

In  addition,  for  every  successful  differential  orbit  correction  a 
punched  card  is  generated,  for  post-processing,  containing  the 
following  information: 


Columns 

Data  Format 

Form 

1-3 

Day  No.  of  the  Year 

13 

XXX 

4-5 

Hour 

12 

XX 

6-7 

Minute 

12 

XX 

8-9 

Seconds 

12 

XX 

10-14 

Geoc.  Lat. 

F5.1 

+XX.X 

15-19 

West  Long. 

F5.1 

XXX. X 

20-26 

Alt.  at  Perigee  (km) 

F7.2 

xxxx.xx 

27-35 

Density  at  Perigee 

E9.3 

X.XXXE-XX 

36-44 

Density  at  Std.  Height 

E9.3 

X. XXXE-XX 

45-51 

Ait.  at  1/2  Scale  Ht.  (km) 

F7.2 

xxxx.xx 

52-60 

Density  at  1/2  Scale  Ht. 

E9.3 

X.XXXE-XX 

61-66 

Model  Factor 

F6.3 

XX. XXX 

67-70 

Time  Span  Used  -  Hours 

F4.1 

XX. X 

71-72 

Atmospheric  Model  Used 

12 

XX 

73-76 

Standard  Error  of 
Differential  Correction 

F4.2 

x.xx 

1.2. 4. 8  BADMEP  Input 


The  run  control  parameters  are  read  through  the  NAMELIST  group  NEWINX 
if  columns  61-64  of  the  device  control  card  contains  a  one,  or  the 
NAMELIST  group  NEWIN  if  it  contains  any  other  value.  The  two  NAMELIST 
groups  contain  the  same  variables  except  that  NEWIN  contains  the 
elements  ELM(1)  -  ELM(6)  and  NEWINX  contains  the  position  and  velocity 
PV0Z(1)  -  PV0Z(6).  The  other  parameters  are  as  described  for  NEWIN  in 
Reference  2,  except  for  the  expansion  of  the  gravity  model  arrays  C  and 
S,  as  in  CADNIP,  and  an  additional  parameter  IDC,  a  numerical  code 
required  by  the  current  AFGL  system  to  indicate  the  plotting  device  to 
be  used.  For  definition  of  this  code  the  user  should  consult  the 
latest  AFGL  User's  Guide  and  systems  bulletins.  The  gravity 

coefficients  Cnm  and  Snm  are  to  be  read  into  the  arrays  C  and  S  as 

follows: 

Cnm:  C  [n(n  +  1 )/2  +  m  -  2j 

Snm:  S  [n(n  -  1)/2  +  m  -  1 J 

The  density  models  available,  and  thpir  identifying  codes,  to  be  input 
to  MODTYP,  are  the  same  as  just  described  for  CADNIP. 
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1.3  Neglect  of  Minor  Constituents 


The  Jacch ia-Bass^,  Jacchia  1977^6,  and  MSIS®>^,31  density  models  are 
detailed  composition  models  requiring  the  computation  of  the  densities 
of  6  or  7  constituents.  If  only  the  total  mass  density  is  desired,  a 
significant  savings  may  be  made  by  omitting  computation  of  minor 
constituents.  A  study  was  made  of  the  relative  importance  of  the 
various  constituents  in  the  Jacchia-Bass  model  for  height  range 
120  km  -  1000  km  and  the  temperature  range  1500  K  -  1900  K  in 
increments  of  10  km  and  100  K.  An  algorithm  was  set  up  to  indicate,  at 
each  height  and  temperature,  which  constituents  could  be  neglected 
leaving  a  total  mass  density  error  of  18.  Based  on  these  results  the 
various  consitutents  may  be  neglected  over  the  following  height  and 
temperature  ranges: 


°2: 

h  >_ 

250  +0.3 

(T-500) 

0  : 

h  >_ 

630  +  1.2 

(T-500) 

N2: 

h  > 

360  +  0.6 

(T-500) 

N  : 

all 

heights 

He: 

h  <_ 

240  +  0.2 

(T-500) 

Ar: 

all 

heights 

H: 

h  <_ 

330  +  0.6 

(T-500) 

where  h  is  the  height  in  km  and  T  is  the  exospheric  temperature  in  K, 
These  modifications  lead  to  a  308  reduction  in  CP  time. 

A  further  138  reduction  was  gained  by  substituting  in-line  code  for 
one-dimensional  interpolation  in  place  of  calls  to  a  general  purpose 
one-and  two-dimensional  interpolation  routine.  Two-dimensional 
interpolation  is  not  used  as  often  as  in  the  original  Oachia  1977  model 
because  tables  are  replaced  by  their  analytic  representation  over  most 
of  the  height  range. 


1 .4  Satellite  Accelerometer  Data  Studies 

Satellite-borne  accelerometera  designed  by  AFGL  scientists  have  been 
the  source  of  plentiful  in-situ  neutral  density  measurments  in  the 
lower  thermosphere.  In  addition  to  recent  new  triaxial  accelerometers, 
flown  during  periods  of  moderate  and  high  solar  activity^,  an  extensive 
data  base  has  been  compiled^  from  four  earlier  designed  single  axis 
systems  flown  on  satellites  S3-1 ,  AE-C,  AE-D,  and  AE-E  during  the 
preceding  solar  cycle.  This  data  base  has  been  used  for  evaluation  of 
recent  density  models  and  development  of  an  empirical  global  density 
model^ .  The  more  recent  triaxial  accelerometers  provide  the 
opportunity  to  extend  the  earlier  studies  to  the  higher  solar  activity 
peak  of  the  current  cycle.  In  contrast  to  the  earlier  instruments, 
they  provide  data  continuously  for  extended  periods  of  time. 
Consequently  more  detailed  studies  can  be  made  of  thermospheric 
response  to  geomagnetic  activity. 

In  addition,  5  months  of  continuous  data  from  a  rotating  single  axis 
accelerometer6  has  became  available  which,  when  supplemented  by  the 
triaxial  accelerometer  data,  makes  possible  detailed  studies  of  the 
latitudinal  behavior  of  the  semiannual  variation  below  220  km. 

The  difficulty  in  precisely  representing  thermospheric  dynamics  stems 
from  a  lack  of  sufficiently  available  and  representative  solar  and 
geomagnetic  indices  as  well  as  a  need  for  understanding  the  physics  of 
this  region.  Studies  emphasizing  the  global  specification  of  densities 
to  meet  specific  accuracy  requirements  and  provide  inputs  for  specific 
missions  are  being  made .  This  entails  statistical  evaluation  of  the 
accuracy  of  existing  models,  using  available  data;  incorporation  of 
simple  mathematical  fixes,  when  appropriate  and  feasible,  to  the 
model(s)  chosen  for  the  density  specification;  and  assessment  of  the 
accuracy  of  the  densities  so  specified. 


Section  1.4.1  will  discuss  software  developed  to  construct  data  bases 
for  use  in  evaluation  of  existing  density  models.  These  include  data, 
the  corresponding  density  model  values,  and  data-to-model  ratios. 
Section  1.4.2  discusses  software  which  uses  these  data  bases  to 
perform  comparative  statistical  evaluations  of  models.  Lastly  Section 
1.4.3  discusses  construction  of  normalized  data  bases  for  study  of  the 
semiannual  variation. 

1.4.1  Data  Bases  for  Density  Model  Evaluation 

AFGL  scientists,  with  analytic  support  from  RMY»  reduce  the  raw 
accelerometer  data  to  densities  and  construct  a  data  base  (Table  4) 
which  includes,  for  each  sample,  the  Greenwich  and  local  times, 
geographic  and  geomagnetic  coordinates,  solar  and  geomagnetic  activity 
indices,  and  3acchia  1971  model  evaluative  information  (model  density 
and  ratio  of  measured  density  to  model  density). 

Program  FOURMOD  has  been  written  to  provide  evaluative  information  for 
up  to  four  CADNIP  models  selected  from  Table  5.  Program  DENDB, 
discussed  later,  performs  a  similar  function  for  the  MSIS  77,  MSIS  78 
and  Jacchia  1977  models.  The  density  model  package  of  CADNIP  has  been 
modified  to  compute  up  to  4  of  these  models  in  parallel  in  a  single 
run.  Due  to  memory  requirements,  a  single  run  can  accommodate  no  more 
than  two  models  which  require  storage  of  the  density  tables  referred  to 
in  Section  1.2. 4. 5  of  this  report  and  Section  A. 3  of  Reference  2.  The 
only  exception  to  this  is  when  models  1  and  2  are  selected;  since  they 
use  the  same  table,  a  third  model  using  a  table  may  then  be  selected. 
The  density  tables  are  assumed  to  be  stored  on  a  master  BCD  file, 
TAPE5,  similar  to  that  used  as  a  "system  file"  by  program  CADNIP 
(Reference  2,  Section  A. 8).  The  user  controls  the  copying  of  selected 
segments  of  this  file  onto  TAPES  1-4  with  a  device  control  card  similar 
to  that  described  in  Reference  2,  Section  A. 8. 


Table  4  Accelerometer  Density  Data  Base 

Header  Record  (1  per  file) 

*  Description 

Word  Count  (40) 

Group  Count  (1) 

Satellite  ID 
Experiment  Name 
Altitude  (km)*  or  Blank** 

Unused 

Data  Records  (1  or  more  per  file) 

Description 

Word  Count  (40) 

Group  Count  (12) 

Orbit  Number 
Date  -  YYDDD 
GMT  -  Total  Seconds 
GMT  -  Hours 
GMT  -  Minutes 
GMT  -  Seconds 
Local  Time  -  Hours 
Local  Time  -  Minutes 
Local  Time  -  Seconds 
Leg  (Ump,  D=down) 

Day/Night  (D  or  N) 

Spun/ Ds pun  (S  or  D) 

Geographic  Latitude  (Deg) 
Geographic  East  Longitude  (Deg) 
Geomagnetic  Latitude  (Deg) 
Geomagnetic  East  Longitude  (Deg) 
Invariant  Latitude  (Deg  )  or  Blank 
Measured  Density  (gm/cm3) 

Jacchia  1971  Model  Density  (gm/cm^) 


Table  4  (Cont'd.) 


,Word 

Description 

Format 

22 

■ 

Ratio  (Meas/Jacchia  1971) 

F 

23,  24 

Unused 

25 

Dally  Average  Ap  or  Blank  *** 

F 

26 

Unused 

27 

Kp  (6.7  hour  lag) 

F 

28 

F10  ?  (1  day  lag) 

F 

29 

Fio  i  (3-solar-rotation  average) 

F 

30 

Altitude  (Km)**  or  Blank* 

F 

31-40 

Unused 

Words  1-40  repeat  for  remaining  11  groups.  Epd-of-file  separates  files; 
double  end-of-file  follows  last  file. 

*  Altitude  Data  Base;  **Time  Data  Base;  ***Da11y  Average  Ap  required  for  computation 

of  MSIS  models  by  program  F0URM0D 


Table  5  Density  Models  Computed  by  Program  F0URM0D 


Indicator 

1 

2 

3 

4 

5 

7 

8 
9 

10 

11 

12 

13 

*  Word  25  of  the  input  database  (Tab! 


Model 

Jacchia  1964 
1966  Supplements 
Jacchia  1971 
USSR  -  Cosmos 
Jacchia-Walker-Bruce 
Lockheed/NASA 
MSIS  77  * 

1962  U.S.  Standard 
MSIS  78  * 

DENSEL 

Jacchia  1970 
Jacchia  1973 

4)  must  contain  the  daily  average  Ap. 


-40- 


This  card  consists  of  8  fields,  read  in  814  format,  specifying  the 
following: 

Field  Description 

1  Number  of  BCD  card  images  to  skip  on 

TAPES  prior  to  copying  to  TAPE1 

2  Number  of  BCD  card  images  to  copy 

subsequently  from  TAPES  to  TAPE1 

3  Number  of  BCD  card  images  to  skip 

subsequently  on  TAPE5,  prior  to 
copying  to  TAPE2 

4  Number  of  BCD  card  images  to  copy 

subsequently  to  TAPE2 

5,6  Similar  to  fields  3  and  4,  but  for 

creating  TAPE 3 

7,8  Similar  to  fields  3  and  4,  but  for 

creating  TAPE4 

Although,  as  previously  noted,  under  current  limitations  at  most  2  of 
the  4  files  so  created  may  subsequently  be  read  by  the  density  model 
package,  FOURMOD  has  provided  4  files  to  permit  later  expansion  if 
memory  restrictions  become  less  severe.  Following  the  device  control 
card,  two  other  cards  are  required  as  described  in  Table  6. 

The  solar  and  geomagnetic  activity  file  normally  used  by  CADNIP  is  not 
required,  since  most  of  the  information  is  already  on  the  input  data 
base  (Table  4).  However  the  six-solar-rotation  average  of  the  daily 
solar  flux  F10.7,  rather  than  the  three-solar -rotation  average,  is  the 
preferred  estimate  of  the  smoothed  solar  flux  for  all  Jacchia  models^, 
in  spite  of  the  reference  to  the  three-solar-rotation  average  in  early 
Jacchia  model  reports  (References  21  and  24).  Therefore  a  separate 
file,  TAPE6,  is  required  containing  one  binary  data  record  consisting 
of  the  following: 


Table  6  Program  FOURMOD  Punched  Card  Input 


Card  # 

Col umn 

Format 

Variable 

Descri pti on 

1 

Device 

Control  Card 

(See  text) 

2 

1-5 

15 

NMOD 

Number  of  models  to  be 
computed  (1,2, 3, or  4) 

2 

6-10 

15 

IM0D(1) 

Indicator  of  first  model 
(see  text) 

2 

11-15 

15 

NDENTP(l) 

Number  of  file  containing 
tables  for  first  model 
(0,1, 2, 3, or  4) 

2 

16-20 

15 

I  MO  D ( 2 ) 

2 

21-25 

15 

NDENTP(2) 

2 

26-30 

15 

IM0D(3) 

2 

31-35 

15 

NDENTP(3) 

2 

36-40 

15 

I MOD ( 4 ) 

2 

41-45 

15 

NDENTP(4) 

3 

1-5 

15 

IREPT 

Repeat  factor  for  printout 
of  data 

3 

6-10 

15 

IALTDB 

Zero  indicates  time  data 
base;  non -zero  indicates 

altitude  data  base 

Zero  indicates  data  words 
3-9,  14,  16  are  integers; 
non-zero  indicates  they  are 
fl oating  point  (real ) 


3 


11-15 


15 


IORMAG 


Word 


Variable 


Description 


1  TSTART  Modified  Julian  Day  of  the 

first  day  for  which  smoothed 
solar  flux  is  provided 

2  N  Length  of  the  period,  in  days, 

for  which  smoothed  solar  flux 
is  provided 

FBAR6  Smoothed  solar  flux  (6-solar- 

rotation  averages  of  daily 
solar  flux)  for  each  day  of 
this  period 

The  data  provided  on  TAPE6  must  of  course  cover  at  least  all  times  for 
which  density  data  are  given  on  the  input  data  base. 

The  input  and  output  data  bases  are  identified  as  TAPE9  and  TAPE10, 
respectively.  Table  7  describes  the  output  data  base. 

The  following  differences  between  the  Jacchia  1971  model  version 
provided  by  FOURMOD  (from  the  CADNIP  package)  and  the  original  model 
(Reference  24)  provided  in  the  input  data  base  should  be  noted. 

1)  In  the  original  model  there  is  a  break  at  200km  in  the 
computation  of  the  geomagnetic  activity  correction:  below 
200km  a  hybrid  correction  is  made  to  the  exospheric 
temperature  and  to  the  density  (Reference  24),  while  above 
200km  only  an  exospheric  temperature  correction  is  computed. 
The  CADNIP  (FOURMOD)  version  uses  the  hybrid  correction  at  all 
heights. 

2)  Seasonal-latitudinal  corrections  for  helium  are  neglected  in 
the  CADNIP  (FOURMOD)  version. 


3  thru 
N  +  2 


Table  7  Program  FOURMOD  Output  Data  Base 


Header  Record  (1  per  file) 

Copy  of  input  header  record 

Data  Records  (1  or  more  per  file) 

Copy  of  input  data  records  with  following  exceptions: 


Words 

Description 

Format 

9,  14,  16 

Same  as  Input 

F 

19 

Density  for  first  model  (gm/cm^) 

F 

20 

Density  for  second  model  (gm/cm^) 

F 

21 

Density  for  third  model  (gm/cm  ) 

F 

22 

Ratio  (Mea s/first  model) 

F 

23 

Ratio  (Meas/second  model) 

F 

24 

Ratio  (Meas/third  model) 

F 

26 

Ap  (6.7  hr  time  lag) 

F 

37 

Density  (fourth  model,  gm/cnr*) 

F 

38 

Ratio  (Meas/fourth  model) 

F 

3)  The  original  371  on  the  input  data  base  is  computed  with  the 
three  solar-rotation  average  of  solar  flux,  while,  as 
mentioned  previously,  FOURMOD  employs  the  six-solar-rotation 
average. 

Of  these  the  first  difference  is  probably  the  most  significant,  since 
the  models  will  differ  above  200km  for  moderate  and  high  geomagnetic 
activity.  The  helium  correction  is  not  expected  to  be  important  in  the 
altitude  range  of  interest.  The  U3e  of  different  solar  flux  averaging 
periods  could  affect  results  for  moderate  to  high  solar  activity. 


Table  8  DENDB  Output  Data  Base 

Header  Record  (1  per  file) 

Copy  of  input  header  record 

Data  Records  (1  or  more  per  file) 

Copy  of  input  data  records,  with  the  following  exceptions: 


Words 

Descri pti on 

Format 

-9,14,16 

Same  as  input 

F 

20 

Jacchia  1977  model  density 

(gm/cirr)  , 

F 

21 

MSIS  77  model  density  (gm/cnr) 

F 

23 

Ratio  (meas/Jacchi a  1977) 

F 

24 

Ratio  (meas/MSIS  77) 

F 

25 

o 

MSIS  78  model  density  (gm/cm  ) 

F 

26 

Ratio  (meas/MSIS  78) 

F 

Program  DENDB  constructs  similar  model  evaluation  data  bases  for  the 
Jacchia  1977,  MSIS  77  and  MSIS  78  density  models.  It  makes  use  of  the 
density  model  computation  package  of  program  DENMOD^.  The  input  data 
base  (Table  4)  must  be  on  TAPE1,  while  the  output  data  base, 
summarized  in  Table  8,  is  written  on  TAPE2.  The  solar  and  geomagnetic 
activity  tables  required  by  CADNIP  are  also  required  by  DENDB  and  must 
be  input  through  TAPE3.  A  mass  storage  file  JSDM,  required  for  the 
Jacchia  1977  model,  is  described  in  Table  1,  Section  1,  of  Reference 
36. 

Density  model  evaluation  data  bases  have  been  constructed  for  the 
models  listed  in  Table  5  and  the  Jacchia  1977  model,  using  data  from 
the  AE/S3-1  data  base,  the  rotatable  calibration  accelerometer  (ROCA) 
and  two  recently-flown  triaxial  accelerometers  (SETA-1  and  SETA-2). 


V 

xd 


DENSITY  MODEL  RATIO  STATISTICAL  SUMMARIES 


NO. 

M1n 

MAX 

MIN 

MAX 

MIN 

1 

O.CCOO 

3.1000 

8. OCOO 

24.0000 

-90.0000 

2 

0.0000 

3. 1000 

8.0000 

24 . OOOO 

-80.0000 

3 

0.0000 

3.1000 

8.0000 

24 .  OOOO 

-70.0000 

4 

0.0000 

3. 1000 

8.0000 

24.0000 

-60.0000 

5 

0.0000 

3. 1000 

8.0000 

24.0000 

-50.0000 

6 

0.0000 

3.1000 

8.0000 

24 . OOOO 

-40.0000 

7 

0.0000 

3. 1000 

8.0000 

24 . OOOO 

-30.0000 

8 

0.0000 

3.1000 

8.0000 

24.0000 

-20.0000 

9 

0.0000 

3.1  COO 

8.0000 

24.0000 

-10.0000 

10 

0.0000 

3.1000 

8.0000 

24 . OOCO 

0.0000 

1  1 

0.0000 

3.  1000 

8.0000 

24 . OOOO 

10. OCOO 

12 

0.0000 

3. 1000 

8.C000 

24 . 000 0 

20.0000 

13 

0.0000 

3.1000 

8.0000 

24 . OOOO 

30.0000 

14 

0.0000 

3. 1000 

8.0000 

24 . OOOO 

40.0000 

ts 

0.0000 

3.1000 

8.0000 

24.0000 

50.0000 

16 

0.0000 

3.1000 

8. OCCO 

24.0000 

60.0000 

17 

0.0000 

3. 1000 

8. OOCO 

24.0000 

70.0000 

IS 

0.0000 

3. 1000 

8. OOCO 

24 . OOCO 

80.0000 

19 

3.1000 

4.4000 

3. OCOO 

24.0000 

-90.000C 

2C 

3. 1000 

4 . 4000 

8.CC00 

24.0000 

-30.0000 

21 

3.1000 

4 . 4000 

8.0000 

24.0000 

-70.0000 

22 

3.1 000 

4 . 4000 

8.0000 

24.0000 

-60.0000 

23 

3. 1000 

4.4000 

8.0000 

24 . OOOO 

-50.0000 

21 

3.1000 

4. 4000 

8.0000 

24.0000 

-40.0000 

25 

3.1000 

4.400C 

8.0000 

24. OOCO 

-30.0000 

20 

3.1000 

4.40 00 

8.0000 

24.0000 

-20.0000 

27 

3.1000 

4.4000 

8-0000 

24.0000 

-10.0000 

23 

3.1000 

4.4000 

8.0000 

2 1 .0000 

0.0000 

23 

3. 10 00 

4.4000 

8. OCCO 

24.0000 

10.0000 

30 

3.1000 

4.4000 

8.0000 

24.0000 

20.0000 

3! 

3.1000 

4.4000 

8. OOCO 

24.0000 

30.0000 

32 

3.  1000 

4.4000 

8.0000 

24 . OOOO 

40.0000 

33 

3.1000 

4.40GC 

8.0000 

24.0000 

50.0000 

34 

3.  1000 

4.4000 

8.0000 

24.0000 

60.000C 

35 

3.1000 

4 .4000 

8.0000 

24.0000 

70.0000 

36 

3.1000 

4.4000 

8.0000 

24.0000 

80.0000 

37 

4.4000 

9.1000 

8.0000 

24.0000 

-90.0000 

33 

4.4000 

9.1000 

8.0000 

24.0000 

-8 0.0000 

34 

4.4000 

9.1000 

8.0000 

24.0000 

-70.0000 

40 

4.4000 

9.1000 

8.COOO 

24.0000 

-60. OOOO 

41 

4.4000 

9.1000 

8.0000 

24.0000 

-50.0000 

42 

4.4000 

9.1000 

8.0000 

24.0000 

-40.0000 

43 

4.4000 

9.1000 

8.0000 

24.0300 

-30.0000 

44 

4. 4000 

9.  1000 

8.0000 

24.0000 

-20.0000 

45 

4.4000 

9.1000 

8.0000 

24.0000 

-10.0000 

40 

4.4000 

9.1000 

8.000C 

24. OOCO 

0.0000 

4’. 

4.4000 

9.  IOCO 

8.0000 

24 . OOOO 

10.0000 

48 

4.4000 

9.1000 

8.0000 

24. OCOO 

20.0000 

43 

4.4000 

9. 1000 

e.ocoo 

24 . OOOO 

30.0000 

50 

4.4000 

9.1000 

B.COOO 

24.0000 

40.0000 

51 

4.4000 

9.1000 

8.0000 

24 .OOCO 

50.0000 

52 

4.4000 

9.1000 

8.0000 

24.0000 

60.0000 

53 

4.4000 

9 . 1  COO 

8.0000 

24.0000 

70.0000 

54 

4.4000 

5.1000 

8.0000 

24.0000 

80.0000 

LAT 

MAX 

-80.0000 
-70.0000 
-60. COCO 
-60.0000 
-40.0000 
-20. COCO 
-20.0000 
-10 .  occo 
o .  oooo 
io.ooco 

20. COCO 
30.0000 
40 . ooco 
50.  OOCO 
60.0000 
70.0000 

go .  oooo 
90.  occo 
-go . oooo 

-70.0000 
-60.0000 
-50.0000 
-4C. OOOO 
-30. OCCO 
-20.0000 
-IP.  0,000 
0 .  POCO 
10.0000 
20.0000 
30.CC00 
40.0000 
50.0000 
CO.  oooo 
70.0000 
60 . OOOO 
90.0000 
-80.0000 
-70. OCCO 
-60.0000 
-50.0COO 
-40.0000 
-30.00.v0 
-20.0090 
-IO.OOC'0 
O.CPOO 

lo.oooo 
20.0000 
30 . coco 
40. OCCO 
50 . ocoo 
60.  O.'CO 
70 .  COCO 
80.0000 
90.0000 


Figure  2.  Program  STAT  Sample:  Three  Dimensional  Bins 


■ 


V.V.V.V.'.-.V. 


Table  9  Punch  Card  Input  for  Program  STAT 


Card  # 

Col umn 

Format 

Varlable(s) 

Descri ptl on 

1 

1 

11 

IALTDB 

0  for  time  data  base 

1  for  altitude  data  base 

'  1 

6-10 

15 

ALTMN 

Minimum  altitude  to  process 

1 

11-15 

15 

AITMX 

Maximum  altitude  to  process 

2 

1-30 

3A10 

Var  (1), 

Var  (2), 

Var  (3) 

BCD  (10  character)  names  of 
up  to  3  variables  to  be  used 
to  define  bins.  If  a  blank 
name  is  given,  only  the  non¬ 
blank  names  preceding  it  wil 
be  considered. 

3 

1-40 

4A10 

OENHED  (1), 

DENHED  (2), 

OENHED  (3), 

DENHED  (4) 

BCD  names  of  the  density 
models  (up  to  4)  to  be 
processed.  A  blank  name 
has  the  same  effect  as  for 
VAR  in  card  2. 

3 

41-60 

415 

IWD<1),  I WD ( 2 ) , 
IWD( 3 ) ,  I W D ( 4 ) 

Word  numbers,  on  data  base 
(Tables  6,8),  where  the 
corresponding  model  ratios 
are  stored. 

4+:  I 

set  of  cards 

per  variable  named  In  card 

2 

4A 

1-5 

15 

IVARN  (I) 

Word  number  on  data  base 
where  the  value  of  the  Ith 
variable  Is  stored 

4A 

6-10 

IS 

N8 IN  (I) 

Number  of  bins  for  this 
variable  (24  maximum) 

4B 

1-80 

8F10.3 

( ( AMINMX ( J  ,  K ,  I ) 
J-1,2)  K=1 ,  NBIN), 
DEFMNMX (1,1), 
DEFMNMX  (2,1) 

AMINMX:  Minimum  and  maximum 
values  of  I th  variable  per 
bin . 

DEFMNMX:  Default  minimum 
and  maximum  values  for  I th 
variable 

5 

1-3 

311 

IBNFLG  (1), 

IBNFLG  (2  , 

IBNFLG  (3) 

Bin  flags  for  the  variables 
named  in  card  2  (see  text). 
Several  cards  may  be  input. 

terminated  by  a  card  of  all 
zeroes . 


Card  46  may  be  continued  as  necessary. 


Table  10.  Sample  Bln  Specification  for  Program  STAT 


Variable 

No.  of  Bins 

Minima 

Maxima 

Default  Minimum 

Default  Maxima 

KP 

3 

0.0 

3.1 

0.0 

9.1 

3.1 

4.4 

4.4 

9.1 

Local  Time 
(Hrs. ) 

1 

8 

24 

8 

24 

Geographic 

Latitude 

18 

-90.0 

_q n  n 

-80.0 

_7n  n 

-90.0 

90.0 

1.4.2  Statistical  Evaluation  of  Density  Models 

Statistical  evaluation  of  density  models  has  been  performed  by 
computation  of  the  mean  M  and  percent  standard  deviation  S  from  the 
mean,  of  the  data-to-model  ratio: 

N 

M  =  E  R./N 
i=l  1 


where 

=  data  to  model  ratio  for  the  ith  sample 
N  =  number  of  samples 

The  data  are  typically  sorted  into  bins  defined  by  1-3  variables  in  the 
data  base. 

In  some  instances,  particularly  when  a  large  set  of  models  is  under 
examination,  only  a  single  overall  evaluation  is  desired. 


The  division  of  data  into  bins  permits  one  to  compare  the  effectiveness 
of  various  models  as  certain  conditions  such  as  latitude,  height,  local 
time  or  geomagnetic  activity  are  varied.  A  significant  variation  of  a 
model's  mean  ratio  over  the  range  of  some  parameter,  such  as  latitude, 
would  indicate  a  deficiency  of  the  model  in  its  dependence  on  that 
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variable.  Large  standard  deviations  often  are  caused  by  failure  to 
include  some  variable  in  the  model.  A  typical  example  is  the  high 
standard  deviations  often  encountered  at  high  latitudes.  This  is 
presumably  caused  by  poor  modeling  of  heating  at  these  latitudes,  or 
what  is  more  often  said,  the  lack  of  a  simple  index,  or  indices,  which 
adequately  represent(s)  the  dynamics  of  such  high  latitude  effects  as 
particle  precipitation  and  joule  heating.  Unmodelled  density  waves  are 
also  sources  of  high  standard  deviations. 

Program  STAT  has  been  developed  to  provide  these  statistical  evaluation 
capabilities.  Table  9  summarizes  the  punched  card  input  required  to 
operate  this  program.  As  indicated  in  the  card  2  input,  up  to  three 
variables  (such  as  geocentric  latitude,  or  local  time)  may  be  used  to 
define  bins.  STAT  will  count  only  those  variables  named  before  the 
first  blank  10  column  field  on  card  2.  If  the  first  field  is  blank 
then  all  the  data  within  the  altitude  limits  defined  on  card  1  is 
handled  as  one  group.  In  this  special  case  there  is  no  binning,  and 
cards  4  and  5  are  not  required.  A  similar  technique  is  used  by  STAT  in 
reading  card  3  to  determine  the  number  of  models  to  be  processed. 
However  in  this  case  a  blank  field  in  columns  1-10  triggers  an  error 
termination  (there  must  be  at  least  one  model). 

A  detailed  explanation  of  input  cards  4  and  5  is  given  here.  The  user 
may  submit  as  many  bin  flag  cards  as  desired,  terminating  with  all 
zeroes.  Each  card  will  generate  all  possible  bins  with  minimum  and 
maximum  pairs  specified  by  the  array  AMINMX  for  each  variable  whose 
flag  is  odd,  and  specified  by  the  array  DEFMNMX  for  each  variable  whose 
flag  is  even.  Suppose,  for  example,  we  are  binning  with  respect  to  KP, 
local  time  and  latitude  and  have  specified  the  bin  minima  and  maxima 
for  each  variable  as  in  Table  10.  Then  the  bin  structure  defined  by 
flags  1,  0,  1  would  be  as  shown  in  Figure  2.  If  the  total  number  of 
bins  constructed  by  all  bin  flag  cards  exceeds  300,  an  error  message  is 
written  and  the  program  terminates. 


The  specification  of  altitude  minimum  and  maximum  on  card  1  appears 
unnecessary  at  this  point.  However  due  to  increasing  uncertainty  in 
the  data  at  high  altitudes  it  is  often  desired  to  limit  the  altitude 
range  of  the  data  to  be  considered,  even  when  the  altitude  is  not  being 
used  as  a  bin  variable.  In  such  cases  the  input  on  card  1  avoids  the 
necessity  of  superficially  using  the  altitude  as  one  of  the  variables 
defining  bins. 

Figure  3  shows  a  sample  printout  of  the  results  for  the  sample  case 
described  in  Table  10.  Some  bins  have  no  samples  and  thus  the  results 
are  zeroed.  The  characteristic  trend  toward  high  standard  deviations 
(high  density  variabiity)  at  high  latitudes  is  evident  for  all  levels 
of  geomagnetic  activity.  Variations  of  the  model  ratios  with  latitude 
and  geomagnetic  activity  are  also  present.  The  Jacchia  1977  model 
exhibits  the  most  pronounced  deviations  in  ratio  from  unity,  indicating 
that  its  response  to  increase  in  geomagnetic  activity  is  higher  than 
the  measured  response  in  both  equatorial  and  polar  regions. 

1.4.3  Studies  of  the  Semiannual  Variation 


The  accelerometer  measurements  present  an  opportunity  for  detailed 
study  of  the  semiannual  variation  at  low  altitudes.  The  Jacchia  1971 
model,  based  mainly  on  analysis  of  satellite  drag  above  200  km, 
predicts  an  altitude  dependence  but  no  latitudinal  dependence.  It 
would  therefore  be  of  interest  to  determine  if  the  lack  of  latitudinal 
dependence  holds  also  for  the  accelerometer  data,  and  if  the  magnitude 
of  the  semiannual  variation  agrees  with  that  extrapolated  downward  from 
the  higher  altitude  results  via  the  Oacchia  1971  model. 
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Figure  3.  Sample  Program  STAT  Output 


The  beat  data  base  constructed  thus  far  is  the  ROCA  data  base,  which 
includes  continous  coverage  from  day  89  to  day  223,  in  1978.  This 
therefore  includes  the  April  maximum  and  the  July  minimum.  Data 
covering  later  portions  of  a  year  may  be  available  shortly. 

Figure  4  summarizes  the  software  system  that  has  been  developed  to 
support  this  study.  Program  BNSORT  accepts  as  input  any  data  base 
prepared  by  Program  FOURMOD  (Table  7)  or  DENDB  (Table  8).  Given  bin 
parameter  and  normalization  specifications  (Table  11),  it  constructs  a 
bin-sorted  data  base  of  densities  (Table  12)  normalized  to 
approximately  remove  all  variations  except  the  semiannual.  The 
normalized  densities  are  computed  by 

^n  =  R  dmod 

where  dn  is  the  normalized  density,  R  is  the  measured  to  model  ratio  on 
the  input  data  base,  and  dmo(i  is  the  model  value  for  the  specified 
normalization  parameters:  altitude,  latitude,  local  time,  solar 
activity  and  geomagnetic  activity,  and  with  the  solar  declination  set 
to  zero  regardless  of  the  time  of  year.  Program  DAILAV  constructs, 
from  this,  an  output  data  base  (Table  13)  containing  the  daily 
averaged  normalized  densities.  Program  SUADB  formats  these  into  a  file 
compatible  with  AFGL's  interactive  graphics  program,  SUATEK. 

Figure  5  shows  some  sample  results.  Despite  the  scatter  in  the  data, 
the  semiannual  variation  appears  to  be  nearly  the  same  at  all 
latitudes.  There  is  some  indication  that  it  may  be  smaller  than  that 
predicted  by  the  Jacchia  1971  model,  as  shown  in  Table  14.  It  must  be 
emphasized  that  these  are  preliminary  findings  subject  to  change  as 
more  data  become  available. 
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2.0  Artifacts  in  SPA  Observations  of  Radio  Wave  Scintillations 


2. 1  Introduction 

The  purpose  of  the  Scintillation  Processor  A  (SPA)  radio  receiving 
system  is  to  investigate  irregularities  in  the  ionosphere  by  examining 
the  properties  of  transionospheric  signals  originating  from  satellites 
in  eccentric  orbits^.  To  provide  a  stable  phase  reference  reguired  for 
extracting  valid  phase  measurements,  an  ultra-stable  local  oscillator 
signal  is  synthesized.  The  frequency  of  this  source  is  shifted  from 
time-to-time  to  accommodate,  within  the  narrow  receiver  bandwidth,  the 
time-varying  Doppler  frequency  shift  characteristic  of  the  received 
signals.  Continuity  of  phase  is  preserved  during  these  frequency 
shifts.  A  programmable  synthesizer  driven  by  a  stable  reference 
oscillator  is  the  hardware  configuration  used  to  generate  these 
time-varying  local  oscillator  injection  signals.  The  basic  receiver 
architecture  is  sketched  in  Figure  1.  The  processing  of  the  received 
signals  is  performed  by  a  software  system  described  in  Reference  1. 

2.2  Contamination  of  SPA  Signals 

In  operation,  a  problem  is  encountered  with  the  SPA  system.  A  'scope 
display  in  the  field  of  the  intensity  of  the  raw  received  signal 
reveals  an  unwanted  sinusoidal  component  riding  atop  the  received 
signal.  Figure  2  illustrates  the  time-domain  manifestation  of  the 
effect.  Here  the  rapidly  varying,  fine-grain  structure  of  the  signal 
intensity  is  caused  by  "beating"  between  the  wanted  and  the 
contaminating  signals.  In  terms  of  the  processed  data,  this  effect  is 
most  clearly  manifested  in  the  signal  phase  and  intensity  spectrograms, 
where  it  appears  as  an  often  prominent  ledge-like  extension  to  the 
roll-off  portion  of  the  spectrum.  Figure  3  illustrates  the  effect  on 
phase  and  intensity  spectra.  An  important  analysis  objective  is  to 
obtain  the  slope  of  the  roll-off  of  the  spectra.  Always  a  nuisance, 
the  contaminant,  when  very  prominent,  can  thwart  the  process  of  slope 
evaluation.  This  unwanted  energy,  dubbed  "coherent  leakage",  is  an 
artifact  which  results  from  the  particular  architecture  employed  in  the 
SPA  receiver. 
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Figure  2.  Effect  of  Coherent  Leakage  in  the  Time  Domain 
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Examination  of  the  schematic  of  the  SPA  receiver  (Figure  1)  reveals  a 
multiple  conversion  superheterodyne  receiver  with  a  nominal  rf  input 
frequency  of  250  MHz  and  a  first  IF  frequency  of  99.16  MHz.  Notice, 
particularly,  that  the  frequencies  50,  250,  and  99.16  MHz  are  all 
generated  (at  relatively  high  levels)  within  the  receiving  system. 
Without  extreme  care  with  respect  to  out-of-band  rejection  for  bandpass 
filters  and  shielding  for  both  components  and  interconnecting  cables, 
there  is  an  invitation  for  these  locally-generated  frequencies  to 
couple  into  the  signal  channel.  The  dashed  lines  in  Figure  1  are 
intended  to  suggest  several  possible  leakage  paths.  Coherent  leakage 
is  the  result  of  some  such  unwanted  coupling.  There  is  some  evidence 
from  the  field  suggesting  that  the  dominant  path  is  one  from  the 
synthesizer  out  to  the  antenna,  coupling  50  MHz  energy  into  the  first 
mixer,  where  a  fifth  harmonic  is  generated.  The  precise  origin, 
though,  is  immaterial:  with  the  basic  mechanism  understood,  the  task 
at  hand  is  its  elimination,  not  by  hardware  modifications,  but  by 
special  signal  processing  operations. 

At  this  point  the  reader  may  wonder  that  leakage  of  an  ultra-stable 
injection  frequency  at  the  nominal  center  frequency  (zero  Hz  reference 
of  the  system)  should  manifest  itself  as  the  frequency  shifted  and 
spread  effect  shown  in  Figure  3.  This  may  be  explained  as  follows: 

Assume  unit  amplitude  for  the  received  signal  (assumed  sinusoidal) 

s(t)  =  sin  wt 

Assume  an  amplitude  a«1  and  frequency  u>+6  for  the  leakage 
i(t)  =  a  sin  {(  w+6  )  t  + 

The  amplitude  of  the  sum  of  these  two  can  be  shown  to  be  approximated 
by 

Amplitude  =  1  +  a  cos  (6t  +  <t>  ) 


i  »■  -  - . 
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A  spectrogram  of  the  amplitude  would  show  a  large  d.c.  component  (the 
wanted  signal)  plus  a  component  at  the  radian  "beat"  frequency,  (the 
leakage).  The  dominant  component  "captures"  the  zero  Hz  position, 
with  the  lesser  component  shifted  in  frequency.  A  spectrogram  of  the 
original  raw  signal  would  reverse  the  situation:  independent  of 
relative  amplitudes,  the  coherent  leakage  component  would  appear  at 
zero  Hz,  and  the  larger  wanted  signal  would  show  the  smearing  effect  of 
a  changing  Doppler  frequency  shift.  It  can  be  understood,  now,  why  the 
spectrally  pure  leakage  appears  spread  and  frequency  shifted:  The 
explanation  lies  in  the  phenomenon  of  capture,  coupled  with  the  Doppler 
drifting  and  periodic  shifting  of  the  signal  frequency.  That  is,  the 
difference  frequency  between  wanted  signal  and  leakage  varies  with 
time;  and  the  wanted  signal  is  concentrated  near  zero  Hz.  Thus  the 
originally  stable  leakage  has  the  smeared  appearance  shown  in  Figure 


A  similar  argument  could  be  applied  to  the  case  of  phase  spectra,  where 
again  the  dominant  component  captures  the  zero  Hz  position. 

2.3  Suppression  of  Leakage 

In  signal  quadrature  component  space,  the  leakage  represents  a 
transformation  of  the  signal  to  a  new  origin.  The  displacement  of  the 
new  origin  from  the  old  is  given  by  the  quadrature  components  of  the 
leakage.  It  will  be  assumed  that  these  components  are  constant  over 
the  signal  analysis  period  (4096  samples  x  1/50  seconds/sample  =  81.92 
secs,  see  Reference  1).  Measurements  indicate  that  the  leakage 
components  are,  in  fact,  reasonably  constant  over  such  a  period. 
However,  there  can  be  appreciable  drift  from  one  analysis  block  to 
another . 

The  second  assumption  to  be  made  is  that  the  received  signal  has  a 
statistical  characterization  such  that  its  temporal  history, 
represented  in  quadrature  space,  has  circular  symmetry  about  the 
origin.  Contributing  to  the  realization  of  this  condition  is  the  fact 
that  the  signal  is  generally  offset  from  the  zero-frequency  reference 
by  some  varying  non-zero  amount,  causing  rotation  of  the  signal's 
phasor  representation  about  the  origin. 


Now,  if  the  statistics  of  the  signal's  phasor  representation  do  not 
vary  with  time,  it  is  intuitively  obvious  that,  as  this  constant 
pattern  undergoes  a  large  number  of  rotations  about  the  origin,  the 
resultant  temporal  history  of  the  signal  will  have  circular  symmetry. 
Therefore,  we  will  make  the  assumption  of  temporal  homogeneity  of  the 
statistics  of  the  signal  over  the  82-second  analysis  period.  This 
condition  appears  to  be  violated  only  rarely;  the  usual  circumstance  is 
one  in  which  the  rate  of  scintillation  is  very  slow  over  the  analysis 
period. 

Figure  4a  suggests  the  appearance  of  the  signal  scatter  plot  in 
quadrature  component  space  under  conditions  of  no  coherent  leakage  and 
no  frequency  offset.  Temporal  homogeneity  of  the  statistics  implies 
that  the  general  pattern  of  this  plot  would  not  change  as  different 
time-blocks  of  data  are  examined.  Figure  4b  suggests  the  effect  of 
having  a  non-zero  frequency  offset:  the  distribution  of  Figure  4a 
rotates  many  times  about  the  origin,  yielding  a  circularly  symmetric 
pattern.  Figure  4c  illustrates  the  effect  of  coherent  leakage:  the 
origin  (center  of  symmetry)  is  shifted.  The  circular  symmetry  of  the 
received  signal  (Figure  4b)  implies  a  zero-mean  process.  The  mean  of 
the  process  shown  in  Figure  4c  is  at  the  center  of  symmetry,  which  is 
displaced  from  the  origin.  The  displacement  represents  the  quadrature 
components  of  the  leakage. 

Thus,  by  forming  the  mean  of  the  observed  signal  (Figure  4c)  over  the 
82-second  data  collections  period,  we  can  form  a  statistical  estimate 
of  the  coherent  leakage  signal.  The  quadrature  components  of  this 
estimate  can  then  be  subtracted  from  the  raw  signal  to  suppress  the 
unwanted  leakage  and  restore  the  received  signal  to  its  natural 
zero-mean  condition. 

This  technique  has  been  implemented  in  the  operational  version  of  the 
APA  processing  software.  Program  TPSCAN  unpacks  all  the  data,  checks 
the  quality  of  data,  and  identifies  blocks  of  data  to  be  processed.  For 
each  such  block,  it  forms  an  estimate  of  the  leakage  component,  which 
is  stored  for  use  in  subsequent  processing  to  suppress  the  leakage. 
Figure  5  compares  spectra  before  and  after  suppression  of  leakage. 
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Oata  analysts  have  an  obligation  to  view  with  suspicion  any  process 
that  involves  modification  of  raw  data  and  which,  as  a  consequence, 
might  affect  the  validity  of  their  deductions. 

Now,  cancellation  of  coherent  leakage  is  accomplished  by  adding  to  each 
of  the  quadrature  components  of  the  raw  signal  a  quantity  that  is 
constant  over  the  analysis  period. 

It  is  of  interest  to  inquire  how  this  procedure  would  affect  the  data 
if  it  were  to  function  incorrectly.  Consider  two  cases:  data  with 
coherent  leakage  present;  and  data  that  is  free  of  leakage. 

Leakage  present.  In  this  case  the  cancellation  procedure,  if  it 
functions  imperfectly,  will  either  not  fully  suppress  the  leakage 
or  it  might  actually  increase  the  apparent  level  of  leakage. 

No  leakage.  In  this  case,  the  algorithm  should  estimate  zero 
leakage  and  should  not  modify  the  data.  If  its  estimate  is  in 
error,  though,  it  will  subtract  this  erroneous  estimate  from  the 
data.  The  result  will  be  data  which  will  have  the  appearance  of 
contamination  by  leakage. 

Thus,  a  malfunction  of  the  cure  is  manifested  as  symptoms  of  the 
disease.  Analysts  should  monitor  processed  data  for  any  such  evidence. 

2. A  References 

1)  Roberts,  F.R.,  "Software  for  Processing  Scintillation  Data  From 
Satellites  in  Eccentric  Orbits,"  AFGL  Technical  Memorandum  No.  81,  April  7, 
1983. 


3.0  Program  PR0PL0K 


3.1  Introduction 


Ionospheric  ducting  and  iono-to-iono  mode  propagation  are  familiar 
phenomena  in  the  field  of  radio  wave  propagation.  Their  occurrence  for 
ground-based  stations,  though,  is  generally  dependent  on  the  existence  of 
special  conditions  such  as  suitably-oriented  electron  density  gradients  in 
the  ionosphere.  Thus,  availability  of  ducted  modes  is  dependent  on  the 
vagaries  of  the  ionosphere.  However,  these  modes  offer  some  hope  of 
affording  coverage  across  regions  throughout  which  the  ionosphere  is 
affected  by  disturbed  conditions,  precluding  use  of  classical  modes.  There 
is  interest,  then,  in  providing  man-made  launching  mechanisms  (such  as 
artificial  field-aligned  ionization)  in  order  to  make  such  modes  available 
when  needed.  Satellite  experiments  are  to  be  performed  to  explore  this 
technology. 

In  the  conduct  of  experiments,  both  classical  and  ducted  modes  will 
generally  be  observed.  Each  will  exhibit  both  a  characteristic  "range" 
delay  time  and  a  Doppler  frequency  shift.  A  method  of  calculating  Doppler 
shift  and  delay  time  will  be  required  for  both  the  planning  of  data 
collection  and  the  analysis  of  measurements.  Program  PR0PL0K  has  been 
developed  to  satisfy  this  requirement.  PR0PL0K  is,  basically,  a  variant  of 
the  satellite  ephemeris  program,  LOKANGL,  described  in  Reference  1. 

Figure  1  illustrates  the  propagation  geometry  for  both  the  ducted  and 
classical  modes.  As  presently  configured,  PR0PL0K  computes  the  range  delay 
for  the  ducted  modes  exclusive  of  the  launching  up-leg  which  illuminates 
field-aligned  ionization  from  which  the  ducted  mode  is  scattered.  For 
purposes  of  calculation,  the  ducted  path  is  assumed  to  extend  from  the 
point  directly  above  the  station  to  the  satellite  location.  The  altitude 
of  the  path  is  assumed  uniform  and  equal  to  the  altitude  of  the  satellite. 
Similarly,  the  reflection  height  for  the  classical  mode(s)  is  assumed  the 
same  for  each  hop  and  is  equal  to  the  satellite  altitude.  PR0P10K  performs 
calculations  for  the  three  lowest  order  modes  capable  of  propagating,  i.e., 
with  elevation  angle  greater  than  some  pre-set  minimum  value  corresponding 
to  the  minimum  radiation  angle  of  the  transmit  antenna  being  used. 


Figure  1.  Propagation  Geometry 


3.2  Functional  Description 


PROPLOK  is  a  modification  of  AFGL's  standard  satellite  ephemeris  program, 
LOKANGL,  tailored  to  provide  HF  propagation  calculations  for  ground 
station-to-satellite  paths,  taking  account  of  the  satellite's  changing 
position  and  velocity  as  it  traverses  its  orbit. 

Major  changes  to  LOKANGL  are  the  following: 

1.  To  reduce  core  requirements  and  permit  use  on  INTERCOM,  several 
unused  subroutines  are  eliminated. 

-  e.g.,  WRSTP,  SOLVIL,  SILL. 

2.  Input  is  totally  revised. 

-  Unchanging  quantities  are  eliminated  from  the  input 
process  and  written  permanently  into  code. 
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-  Element  sets  are  input  by  means  of  a  chronologically 
ordered  permanent  file.  The  latter  is  to  be  updated 
using  NETED  as  new  element  sets  are  received. 

-  User  cues  added. 

-  Essentially,  only  analysis  time-window  is  required  as 
input  (start  date,  time;  stop  date,  time). 

3.  Propagation  calculations  are  inserted  in  Subroutine  SPPROU. 

4.  Almost  all  standard  LOKANGL  output  is  suppressed,  and  new  output 
is  provided. 

The  basic  LOKANGL  package  is  used  to  furnish  three  types  of  data  for 
propagation  analysis: 

1.  Satellite  coordinates  (as  function  of  time). 

-  These,  together  with  ground  station  locations,  specify 
the  propagation  path  geometry. 

2.  Satellite  velocity  components  in  station-centered  spherical 
coordinate  system  (range,  azimuth,  elevation). 

-  With  vehicle  velocity  and  ray  path  orientation  at  the 
satellite's  location  both  expressed  in  the  same 
coordinate  frame,  the  Doppler  shift  of  signals  can  be 
calculated. 

3.  Bearing  of  satellite  from  ground  station. 

In  providing  the  foregoing  ephemeris  data,  PROPLOK  functions  in  the  normal 
LOKANGL  fashion,  as  described  in  Reference  1.  Next,  the  ephemeris  data  is 
fed  to  the  propagation  calculations,  the  results  of  which  are  printed  out. 
These  functions  are  performed  within  subroutine  SPPROU. 

Figure  2  is  a  modified  version  of  Figure  1  of  Reference  1.  Functions  not 
performed  in  PROPLOK  have  been  eliminated;  new  PROPLOK  functions  have  been 
added  to  the  drawing. 


Read 

Pos/Vel 

Vectors 


Figure  2.  Simplified  Operational  and  Data  Flow  for  PROPLOK 


The  user  is  required  to  furnish  minimal  input  data  to  PROPLOK  consisting 
of: 


Time  interval  (seconds)  between  successive  calculations 
Start  time  of  run  (year,  month,  day,  hour,  minute,  second) 
End  of  time  of  run  (year,  month,  day,  hour,  minute,  second) 


These  data  are  inputted  interactively  in  accordance  with  cues  which  are 
furnished. 


Element  sets  for  the  satellite  are  obtained  from  a  chronologically  ordered 
file  which  must  be  attached  for  each  run.  PROPLOK  scans  this  file  and 
selects  the  appropriate  elements. 

3.3  Input  and  Output 

The  user  is  required  to  furnish  minimal  input  data  to  PROPLOK  consisting 
of : 

-  Time  interval  (seconds)  between  successive  calculations 

-  Start  time  of  run  (year,  month,  day,  hour,  minute,  second) 

-  End  time  of  run  (year,  month,  day,  hour,  minute,  second) 


These  data  are  inputted  interactively  in  accordance  with  cues  which  are 
furnished. 

Element  sets  for  the  satellite  are  obtained  from  a  chronologically  ordered 
file,  TAPE8  (ELSETFILE),  which  must  be  attached  for  each  run.  PROPLOK 
scans  this  file  and  selects  the  appropriate  elements. 

It  is  assumed  that  orbital  data  for  the  satellite  which  carries  the 
instrumentation  for  the  ducted  mode  experiments  will  be  furnished  as  SCF 
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position/velocity  vectors,  records  in  the  TAPE8  file  consist  of  card 
images  of  these  SCF  cards,  the  format  of  which  is  given  on  page  62  of 
Reference  1 . 


An  example  of  the  output  from  PROPLOK  is  presented  in  Figure  3. 
3.4  Mathematical  Approach 

3.4.1  Propagation  Geometries 

3. 4. 1.1  Classical  Modes 


The  propagation  geometry  depends  upon  the  coordinates  of  the  ground  station 
and  the  satellite.  The  former  are  fixed  and  are  imbedded  in  the  code.  The 
latter  are  calculated  in  PROPLOK  as  described  in  Reference  1. 

If 


D  =  Total  ground  distance  of  path 


and 


C  =  Central  angle  subtended  by  total  path 


then 


C  =  D/R  =  arc  cos  (sin  (LAT1)  sin  (LAT2)  + 

cos  (LAT1 )  cos  (LAT2)  cos  (L0N2  -  L0N1))  (1) 

where  R  =  Radius  of  earth 

LAT1,  L0N1  =  Latitude  and  longitude  of  station 

LAT2,  L0N2  =  Latitude  and  longitude  of  satellite. 
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Figure  3.  Sample  Output  From  PROPLOK 


Since  arc  cos  returns  an  angle  in  the  range  (0,  180°),  we  are  assured  that 
the  results  apply  to  the  "short"  path.  We  assume: 

-  the  propagation  path  is  comprised  of  N  half-hops,  where  N  is 
always  an  odd  integer. 

-  ionospheric  reflection  heights  along  the  path  and  satellite 
altitude  are  all  equal  and  are  denoted  H. 

Examining  the  geometry  of  Figure  4,  we  observe  that,  for  right-triangle 
QAB: 

P  =  R  sin  (C/N)  /  cos  (E  +  C/N).  (2) 

Here 

P  =  Propagation  path  length  of  single  half-hop 
C  =  Central  angle  of  total  path 
C/N  =  Central  angle  of  a  half-hop 
R  =  Radius  of  earth 
H  =  Height  of  ionosphere. 

Referring  again  to  Figure  4,  we  can  also  write 

P  sin(E  +  C/N)  =  H  +  R(1  -  cos(C/N))  (3) 

where  E  is  the  elevation  angle  of  the  radio  wave.  Eliminating  P  from  (2) 
and  (3),  we  have: 

E  =  arc  tan((H  +  R(1  -  cos(C/N)))  /  R  sin(C/N))  -  C/N.  (4) 

Knowing  the  coordinates  of  both  the  station  and  the  satellite,  we  can  use 
these  relations  to  solve  for  the  key  propagation  parameters,  E  and  P,  for 
the  classical  modes. 

The  PR0PL0K  calculation  begins  by  assuming  N  =  1  (where  N  is  the  number  of 
half  hops)  and  calculates  E.  N  is  then  incremented  by  2  and  the  solution 
repeated.  The  process  repeats  until  E  is  equal  to  or  greater  than  a 
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specified  minimum  value,  which  typically  could  be  zero  or  the  minimum 
effective  take-off  angle  for  the  particular  antenna  used  at  the  station. 
We  have  now  identified  the  lowest  order  propagating  mode.  N  is  then 
incremented  twice  more  to  provide  calculations  for  each  of  the  three  lowest 
order  propagating  modes. 

3. 4. 1.2  Ducted  Mode 


The  propagation  path  for  the  ducted  mode  is  taken  to  be  a  circular  arc, 
concentric  with  the  earth,  at  the  uniform  ionospheric  height  H.  This  is 
illustrated  in  Figure  5.  This  arc  extends  from  the  point  directly  above 
the  station  to  the  satellite.  The  propagation  path  is  given  by: 

PD  =  (R  +  H)  •  C 

3.4.2.  Propagation  Range  Rate 

The  Doppler  frequency  shift  imparted  by  motion  of  the  satellite  is  obtained 
from  the  time  derivative  of  propagation  range.  This  quantity  is  determined 
from  the  scalar  product  of  vehicle  velocity  and  the  unit  vector  at  the 
satellite  location  aligned  in  the  direction  of  propagation.  LOKANGL 
provides  vehicle  velocity  expressed  in  terms  of  station  look  angle 
coordinates:  time  rates  of  change  of  (station)  range,  azimuth,  and 

elevation.  The  azimuthal  component  of  vehicle  velocity  is  transverse  to 
the  plane  of  propagation  and,  therefore,  makes  no  contribution.  The 
propagation  range  rate  is  found,  then,  by  combining  the  contributions  of 
the  range  and  elevation  components. 

3.4. 2.1  Classical  Modes 

The  geometry  for  propagation  range  rate  calculation  is  shown  in  Figure  6. 
The  calculation  requires  the  value  of  angle  B.  Observe  that  angle  D  is 
given  by: 


D  =  90°  -  C  -  EL 


where  EL  is  the  elevation  look  angle.  Now 
B  =  0  -  (90  -  G)  . 

From  Figure  4  showing  the  half-hop  geometry,  it  may  be  observed  that 
G  =  E  +  C/N  . 

Therefore 

B  s  E  -  EL  -  C  (1  -  1/N) 

The  radial  component  of  velocity  makes  the  contribution 
VR  cos  B  , 

and  the  elevation  contribution  is 

VT  sin  B  . 

3. 4. 2. 2  Ducted  Mode 

The  geometry  for  this  case  is  shown  in  Figure  7.  The  radisl  component 
makes  the  contribution 

VR  cos  (90°  -  G)  =  VR  sin  G  , 

and  the  elevation  contribution  is 

-  VT  cos  G 

The  range  rate  is  given  by  their  sum. 

3.5  References 

1)  Bass,  3.N.,  et  al,  "Analysis  and  Programming  for  Research  in  the 
Physics  of  the  Upper  Atmosphere,"  AFGL-TR-81 -0293,  Logicon  Final 
Report,  October  9,  1981. 


